Line data Source code
1 : /*
2 : * Copyright 2008-2009 Katholieke Universiteit Leuven
3 : * Copyright 2010 INRIA Saclay
4 : * Copyright 2016-2017 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 : #include <isl_ctx_private.h>
15 : #include "isl_map_private.h"
16 : #include <isl_seq.h>
17 : #include "isl_tab.h"
18 : #include "isl_sample.h"
19 : #include <isl_mat_private.h>
20 : #include <isl_vec_private.h>
21 : #include <isl_aff_private.h>
22 : #include <isl_constraint_private.h>
23 : #include <isl_options_private.h>
24 : #include <isl_config.h>
25 :
26 : #include <bset_to_bmap.c>
27 :
28 : /*
29 : * The implementation of parametric integer linear programming in this file
30 : * was inspired by the paper "Parametric Integer Programming" and the
31 : * report "Solving systems of affine (in)equalities" by Paul Feautrier
32 : * (and others).
33 : *
34 : * The strategy used for obtaining a feasible solution is different
35 : * from the one used in isl_tab.c. In particular, in isl_tab.c,
36 : * upon finding a constraint that is not yet satisfied, we pivot
37 : * in a row that increases the constant term of the row holding the
38 : * constraint, making sure the sample solution remains feasible
39 : * for all the constraints it already satisfied.
40 : * Here, we always pivot in the row holding the constraint,
41 : * choosing a column that induces the lexicographically smallest
42 : * increment to the sample solution.
43 : *
44 : * By starting out from a sample value that is lexicographically
45 : * smaller than any integer point in the problem space, the first
46 : * feasible integer sample point we find will also be the lexicographically
47 : * smallest. If all variables can be assumed to be non-negative,
48 : * then the initial sample value may be chosen equal to zero.
49 : * However, we will not make this assumption. Instead, we apply
50 : * the "big parameter" trick. Any variable x is then not directly
51 : * used in the tableau, but instead it is represented by another
52 : * variable x' = M + x, where M is an arbitrarily large (positive)
53 : * value. x' is therefore always non-negative, whatever the value of x.
54 : * Taking as initial sample value x' = 0 corresponds to x = -M,
55 : * which is always smaller than any possible value of x.
56 : *
57 : * The big parameter trick is used in the main tableau and
58 : * also in the context tableau if isl_context_lex is used.
59 : * In this case, each tableaus has its own big parameter.
60 : * Before doing any real work, we check if all the parameters
61 : * happen to be non-negative. If so, we drop the column corresponding
62 : * to M from the initial context tableau.
63 : * If isl_context_gbr is used, then the big parameter trick is only
64 : * used in the main tableau.
65 : */
66 :
67 : struct isl_context;
68 : struct isl_context_op {
69 : /* detect nonnegative parameters in context and mark them in tab */
70 : struct isl_tab *(*detect_nonnegative_parameters)(
71 : struct isl_context *context, struct isl_tab *tab);
72 : /* return temporary reference to basic set representation of context */
73 : struct isl_basic_set *(*peek_basic_set)(struct isl_context *context);
74 : /* return temporary reference to tableau representation of context */
75 : struct isl_tab *(*peek_tab)(struct isl_context *context);
76 : /* add equality; check is 1 if eq may not be valid;
77 : * update is 1 if we may want to call ineq_sign on context later.
78 : */
79 : void (*add_eq)(struct isl_context *context, isl_int *eq,
80 : int check, int update);
81 : /* add inequality; check is 1 if ineq may not be valid;
82 : * update is 1 if we may want to call ineq_sign on context later.
83 : */
84 : void (*add_ineq)(struct isl_context *context, isl_int *ineq,
85 : int check, int update);
86 : /* check sign of ineq based on previous information.
87 : * strict is 1 if saturation should be treated as a positive sign.
88 : */
89 : enum isl_tab_row_sign (*ineq_sign)(struct isl_context *context,
90 : isl_int *ineq, int strict);
91 : /* check if inequality maintains feasibility */
92 : int (*test_ineq)(struct isl_context *context, isl_int *ineq);
93 : /* return index of a div that corresponds to "div" */
94 : int (*get_div)(struct isl_context *context, struct isl_tab *tab,
95 : struct isl_vec *div);
96 : /* insert div "div" to context at "pos" and return non-negativity */
97 : isl_bool (*insert_div)(struct isl_context *context, int pos,
98 : __isl_keep isl_vec *div);
99 : int (*detect_equalities)(struct isl_context *context,
100 : struct isl_tab *tab);
101 : /* return row index of "best" split */
102 : int (*best_split)(struct isl_context *context, struct isl_tab *tab);
103 : /* check if context has already been determined to be empty */
104 : int (*is_empty)(struct isl_context *context);
105 : /* check if context is still usable */
106 : int (*is_ok)(struct isl_context *context);
107 : /* save a copy/snapshot of context */
108 : void *(*save)(struct isl_context *context);
109 : /* restore saved context */
110 : void (*restore)(struct isl_context *context, void *);
111 : /* discard saved context */
112 : void (*discard)(void *);
113 : /* invalidate context */
114 : void (*invalidate)(struct isl_context *context);
115 : /* free context */
116 : __isl_null struct isl_context *(*free)(struct isl_context *context);
117 : };
118 :
119 : /* Shared parts of context representation.
120 : *
121 : * "n_unknown" is the number of final unknown integer divisions
122 : * in the input domain.
123 : */
124 : struct isl_context {
125 : struct isl_context_op *op;
126 : int n_unknown;
127 : };
128 :
129 : struct isl_context_lex {
130 : struct isl_context context;
131 : struct isl_tab *tab;
132 : };
133 :
134 : /* A stack (linked list) of solutions of subtrees of the search space.
135 : *
136 : * "ma" describes the solution as a function of "dom".
137 : * In particular, the domain space of "ma" is equal to the space of "dom".
138 : *
139 : * If "ma" is NULL, then there is no solution on "dom".
140 : */
141 : struct isl_partial_sol {
142 : int level;
143 : struct isl_basic_set *dom;
144 : isl_multi_aff *ma;
145 :
146 : struct isl_partial_sol *next;
147 : };
148 :
149 : struct isl_sol;
150 : struct isl_sol_callback {
151 : struct isl_tab_callback callback;
152 : struct isl_sol *sol;
153 : };
154 :
155 : /* isl_sol is an interface for constructing a solution to
156 : * a parametric integer linear programming problem.
157 : * Every time the algorithm reaches a state where a solution
158 : * can be read off from the tableau, the function "add" is called
159 : * on the isl_sol passed to find_solutions_main. In a state where
160 : * the tableau is empty, "add_empty" is called instead.
161 : * "free" is called to free the implementation specific fields, if any.
162 : *
163 : * "error" is set if some error has occurred. This flag invalidates
164 : * the remainder of the data structure.
165 : * If "rational" is set, then a rational optimization is being performed.
166 : * "level" is the current level in the tree with nodes for each
167 : * split in the context.
168 : * If "max" is set, then a maximization problem is being solved, rather than
169 : * a minimization problem, which means that the variables in the
170 : * tableau have value "M - x" rather than "M + x".
171 : * "n_out" is the number of output dimensions in the input.
172 : * "space" is the space in which the solution (and also the input) lives.
173 : *
174 : * The context tableau is owned by isl_sol and is updated incrementally.
175 : *
176 : * There are currently two implementations of this interface,
177 : * isl_sol_map, which simply collects the solutions in an isl_map
178 : * and (optionally) the parts of the context where there is no solution
179 : * in an isl_set, and
180 : * isl_sol_pma, which collects an isl_pw_multi_aff instead.
181 : */
182 : struct isl_sol {
183 : int error;
184 : int rational;
185 : int level;
186 : int max;
187 : int n_out;
188 : isl_space *space;
189 : struct isl_context *context;
190 : struct isl_partial_sol *partial;
191 : void (*add)(struct isl_sol *sol,
192 : __isl_take isl_basic_set *dom, __isl_take isl_multi_aff *ma);
193 : void (*add_empty)(struct isl_sol *sol, struct isl_basic_set *bset);
194 : void (*free)(struct isl_sol *sol);
195 : struct isl_sol_callback dec_level;
196 : };
197 :
198 0 : static void sol_free(struct isl_sol *sol)
199 : {
200 : struct isl_partial_sol *partial, *next;
201 0 : if (!sol)
202 0 : return;
203 0 : for (partial = sol->partial; partial; partial = next) {
204 0 : next = partial->next;
205 0 : isl_basic_set_free(partial->dom);
206 0 : isl_multi_aff_free(partial->ma);
207 0 : free(partial);
208 : }
209 0 : isl_space_free(sol->space);
210 0 : if (sol->context)
211 0 : sol->context->op->free(sol->context);
212 0 : sol->free(sol);
213 0 : free(sol);
214 : }
215 :
216 : /* Push a partial solution represented by a domain and function "ma"
217 : * onto the stack of partial solutions.
218 : * If "ma" is NULL, then "dom" represents a part of the domain
219 : * with no solution.
220 : */
221 0 : static void sol_push_sol(struct isl_sol *sol,
222 : __isl_take isl_basic_set *dom, __isl_take isl_multi_aff *ma)
223 : {
224 : struct isl_partial_sol *partial;
225 :
226 0 : if (sol->error || !dom)
227 : goto error;
228 :
229 0 : partial = isl_alloc_type(dom->ctx, struct isl_partial_sol);
230 0 : if (!partial)
231 0 : goto error;
232 :
233 0 : partial->level = sol->level;
234 0 : partial->dom = dom;
235 0 : partial->ma = ma;
236 0 : partial->next = sol->partial;
237 :
238 0 : sol->partial = partial;
239 :
240 0 : return;
241 : error:
242 0 : isl_basic_set_free(dom);
243 0 : isl_multi_aff_free(ma);
244 0 : sol->error = 1;
245 : }
246 :
247 : /* Check that the final columns of "M", starting at "first", are zero.
248 : */
249 0 : static isl_stat check_final_columns_are_zero(__isl_keep isl_mat *M,
250 : unsigned first)
251 : {
252 : int i;
253 : unsigned rows, cols, n;
254 :
255 0 : if (!M)
256 0 : return isl_stat_error;
257 0 : rows = isl_mat_rows(M);
258 0 : cols = isl_mat_cols(M);
259 0 : n = cols - first;
260 0 : for (i = 0; i < rows; ++i)
261 0 : if (isl_seq_first_non_zero(M->row[i] + first, n) != -1)
262 0 : isl_die(isl_mat_get_ctx(M), isl_error_internal,
263 : "final columns should be zero",
264 : return isl_stat_error);
265 0 : return isl_stat_ok;
266 : }
267 :
268 : /* Set the affine expressions in "ma" according to the rows in "M", which
269 : * are defined over the local space "ls".
270 : * The matrix "M" may have extra (zero) columns beyond the number
271 : * of variables in "ls".
272 : */
273 0 : static __isl_give isl_multi_aff *set_from_affine_matrix(
274 : __isl_take isl_multi_aff *ma, __isl_take isl_local_space *ls,
275 : __isl_take isl_mat *M)
276 : {
277 : int i, dim;
278 : isl_aff *aff;
279 :
280 0 : if (!ma || !ls || !M)
281 : goto error;
282 :
283 0 : dim = isl_local_space_dim(ls, isl_dim_all);
284 0 : if (check_final_columns_are_zero(M, 1 + dim) < 0)
285 0 : goto error;
286 0 : for (i = 1; i < M->n_row; ++i) {
287 0 : aff = isl_aff_alloc(isl_local_space_copy(ls));
288 0 : if (aff) {
289 0 : isl_int_set(aff->v->el[0], M->row[0][0]);
290 0 : isl_seq_cpy(aff->v->el + 1, M->row[i], 1 + dim);
291 : }
292 0 : aff = isl_aff_normalize(aff);
293 0 : ma = isl_multi_aff_set_aff(ma, i - 1, aff);
294 : }
295 0 : isl_local_space_free(ls);
296 0 : isl_mat_free(M);
297 :
298 0 : return ma;
299 : error:
300 0 : isl_local_space_free(ls);
301 0 : isl_mat_free(M);
302 0 : isl_multi_aff_free(ma);
303 0 : return NULL;
304 : }
305 :
306 : /* Push a partial solution represented by a domain and mapping M
307 : * onto the stack of partial solutions.
308 : *
309 : * The affine matrix "M" maps the dimensions of the context
310 : * to the output variables. Convert it into an isl_multi_aff and
311 : * then call sol_push_sol.
312 : *
313 : * Note that the description of the initial context may have involved
314 : * existentially quantified variables, in which case they also appear
315 : * in "dom". These need to be removed before creating the affine
316 : * expression because an affine expression cannot be defined in terms
317 : * of existentially quantified variables without a known representation.
318 : * Since newly added integer divisions are inserted before these
319 : * existentially quantified variables, they are still in the final
320 : * positions and the corresponding final columns of "M" are zero
321 : * because align_context_divs adds the existentially quantified
322 : * variables of the context to the main tableau without any constraints and
323 : * any equality constraints that are added later on can only serve
324 : * to eliminate these existentially quantified variables.
325 : */
326 0 : static void sol_push_sol_mat(struct isl_sol *sol,
327 : __isl_take isl_basic_set *dom, __isl_take isl_mat *M)
328 : {
329 : isl_local_space *ls;
330 : isl_multi_aff *ma;
331 : int n_div, n_known;
332 :
333 0 : n_div = isl_basic_set_dim(dom, isl_dim_div);
334 0 : n_known = n_div - sol->context->n_unknown;
335 :
336 0 : ma = isl_multi_aff_alloc(isl_space_copy(sol->space));
337 0 : ls = isl_basic_set_get_local_space(dom);
338 0 : ls = isl_local_space_drop_dims(ls, isl_dim_div,
339 0 : n_known, n_div - n_known);
340 0 : ma = set_from_affine_matrix(ma, ls, M);
341 :
342 0 : if (!ma)
343 0 : dom = isl_basic_set_free(dom);
344 0 : sol_push_sol(sol, dom, ma);
345 0 : }
346 :
347 : /* Pop one partial solution from the partial solution stack and
348 : * pass it on to sol->add or sol->add_empty.
349 : */
350 0 : static void sol_pop_one(struct isl_sol *sol)
351 : {
352 : struct isl_partial_sol *partial;
353 :
354 0 : partial = sol->partial;
355 0 : sol->partial = partial->next;
356 :
357 0 : if (partial->ma)
358 0 : sol->add(sol, partial->dom, partial->ma);
359 : else
360 0 : sol->add_empty(sol, partial->dom);
361 0 : free(partial);
362 0 : }
363 :
364 : /* Return a fresh copy of the domain represented by the context tableau.
365 : */
366 0 : static struct isl_basic_set *sol_domain(struct isl_sol *sol)
367 : {
368 : struct isl_basic_set *bset;
369 :
370 0 : if (sol->error)
371 0 : return NULL;
372 :
373 0 : bset = isl_basic_set_dup(sol->context->op->peek_basic_set(sol->context));
374 0 : bset = isl_basic_set_update_from_tab(bset,
375 0 : sol->context->op->peek_tab(sol->context));
376 :
377 0 : return bset;
378 : }
379 :
380 : /* Check whether two partial solutions have the same affine expressions.
381 : */
382 0 : static isl_bool same_solution(struct isl_partial_sol *s1,
383 : struct isl_partial_sol *s2)
384 : {
385 0 : if (!s1->ma != !s2->ma)
386 0 : return isl_bool_false;
387 0 : if (!s1->ma)
388 0 : return isl_bool_true;
389 :
390 0 : return isl_multi_aff_plain_is_equal(s1->ma, s2->ma);
391 : }
392 :
393 : /* Swap the initial two partial solutions in "sol".
394 : *
395 : * That is, go from
396 : *
397 : * sol->partial = p1; p1->next = p2; p2->next = p3
398 : *
399 : * to
400 : *
401 : * sol->partial = p2; p2->next = p1; p1->next = p3
402 : */
403 0 : static void swap_initial(struct isl_sol *sol)
404 : {
405 : struct isl_partial_sol *partial;
406 :
407 0 : partial = sol->partial;
408 0 : sol->partial = partial->next;
409 0 : partial->next = partial->next->next;
410 0 : sol->partial->next = partial;
411 0 : }
412 :
413 : /* Combine the initial two partial solution of "sol" into
414 : * a partial solution with the current context domain of "sol" and
415 : * the function description of the second partial solution in the list.
416 : * The level of the new partial solution is set to the current level.
417 : *
418 : * That is, the first two partial solutions (D1,M1) and (D2,M2) are
419 : * replaced by (D,M2), where D is the domain of "sol", which is assumed
420 : * to be the union of D1 and D2, while M1 is assumed to be equal to M2
421 : * (at least on D1).
422 : */
423 0 : static isl_stat combine_initial_into_second(struct isl_sol *sol)
424 : {
425 : struct isl_partial_sol *partial;
426 : isl_basic_set *bset;
427 :
428 0 : partial = sol->partial;
429 :
430 0 : bset = sol_domain(sol);
431 0 : isl_basic_set_free(partial->next->dom);
432 0 : partial->next->dom = bset;
433 0 : partial->next->level = sol->level;
434 :
435 0 : if (!bset)
436 0 : return isl_stat_error;
437 :
438 0 : sol->partial = partial->next;
439 0 : isl_basic_set_free(partial->dom);
440 0 : isl_multi_aff_free(partial->ma);
441 0 : free(partial);
442 :
443 0 : return isl_stat_ok;
444 : }
445 :
446 : /* Are "ma1" and "ma2" equal to each other on "dom"?
447 : *
448 : * Combine "ma1" and "ma2" with "dom" and check if the results are the same.
449 : * "dom" may have existentially quantified variables. Eliminate them first
450 : * as otherwise they would have to be eliminated twice, in a more complicated
451 : * context.
452 : */
453 0 : static isl_bool equal_on_domain(__isl_keep isl_multi_aff *ma1,
454 : __isl_keep isl_multi_aff *ma2, __isl_keep isl_basic_set *dom)
455 : {
456 : isl_set *set;
457 : isl_pw_multi_aff *pma1, *pma2;
458 : isl_bool equal;
459 :
460 0 : set = isl_basic_set_compute_divs(isl_basic_set_copy(dom));
461 0 : pma1 = isl_pw_multi_aff_alloc(isl_set_copy(set),
462 : isl_multi_aff_copy(ma1));
463 0 : pma2 = isl_pw_multi_aff_alloc(set, isl_multi_aff_copy(ma2));
464 0 : equal = isl_pw_multi_aff_is_equal(pma1, pma2);
465 0 : isl_pw_multi_aff_free(pma1);
466 0 : isl_pw_multi_aff_free(pma2);
467 :
468 0 : return equal;
469 : }
470 :
471 : /* The initial two partial solutions of "sol" are known to be at
472 : * the same level.
473 : * If they represent the same solution (on different parts of the domain),
474 : * then combine them into a single solution at the current level.
475 : * Otherwise, pop them both.
476 : *
477 : * Even if the two partial solution are not obviously the same,
478 : * one may still be a simplification of the other over its own domain.
479 : * Also check if the two sets of affine functions are equal when
480 : * restricted to one of the domains. If so, combine the two
481 : * using the set of affine functions on the other domain.
482 : * That is, for two partial solutions (D1,M1) and (D2,M2),
483 : * if M1 = M2 on D1, then the pair of partial solutions can
484 : * be replaced by (D1+D2,M2) and similarly when M1 = M2 on D2.
485 : */
486 0 : static isl_stat combine_initial_if_equal(struct isl_sol *sol)
487 : {
488 : struct isl_partial_sol *partial;
489 : isl_bool same;
490 :
491 0 : partial = sol->partial;
492 :
493 0 : same = same_solution(partial, partial->next);
494 0 : if (same < 0)
495 0 : return isl_stat_error;
496 0 : if (same)
497 0 : return combine_initial_into_second(sol);
498 0 : if (partial->ma && partial->next->ma) {
499 0 : same = equal_on_domain(partial->ma, partial->next->ma,
500 0 : partial->dom);
501 0 : if (same < 0)
502 0 : return isl_stat_error;
503 0 : if (same)
504 0 : return combine_initial_into_second(sol);
505 0 : same = equal_on_domain(partial->ma, partial->next->ma,
506 0 : partial->next->dom);
507 0 : if (same) {
508 0 : swap_initial(sol);
509 0 : return combine_initial_into_second(sol);
510 : }
511 : }
512 :
513 0 : sol_pop_one(sol);
514 0 : sol_pop_one(sol);
515 :
516 0 : return isl_stat_ok;
517 : }
518 :
519 : /* Pop all solutions from the partial solution stack that were pushed onto
520 : * the stack at levels that are deeper than the current level.
521 : * If the two topmost elements on the stack have the same level
522 : * and represent the same solution, then their domains are combined.
523 : * This combined domain is the same as the current context domain
524 : * as sol_pop is called each time we move back to a higher level.
525 : * If the outer level (0) has been reached, then all partial solutions
526 : * at the current level are also popped off.
527 : */
528 0 : static void sol_pop(struct isl_sol *sol)
529 : {
530 : struct isl_partial_sol *partial;
531 :
532 0 : if (sol->error)
533 0 : return;
534 :
535 0 : partial = sol->partial;
536 0 : if (!partial)
537 0 : return;
538 :
539 0 : if (partial->level == 0 && sol->level == 0) {
540 0 : for (partial = sol->partial; partial; partial = sol->partial)
541 0 : sol_pop_one(sol);
542 0 : return;
543 : }
544 :
545 0 : if (partial->level <= sol->level)
546 0 : return;
547 :
548 0 : if (partial->next && partial->next->level == partial->level) {
549 0 : if (combine_initial_if_equal(sol) < 0)
550 0 : goto error;
551 : } else
552 0 : sol_pop_one(sol);
553 :
554 0 : if (sol->level == 0) {
555 0 : for (partial = sol->partial; partial; partial = sol->partial)
556 0 : sol_pop_one(sol);
557 0 : return;
558 : }
559 :
560 : if (0)
561 0 : error: sol->error = 1;
562 : }
563 :
564 0 : static void sol_dec_level(struct isl_sol *sol)
565 : {
566 0 : if (sol->error)
567 0 : return;
568 :
569 0 : sol->level--;
570 :
571 0 : sol_pop(sol);
572 : }
573 :
574 0 : static isl_stat sol_dec_level_wrap(struct isl_tab_callback *cb)
575 : {
576 0 : struct isl_sol_callback *callback = (struct isl_sol_callback *)cb;
577 :
578 0 : sol_dec_level(callback->sol);
579 :
580 0 : return callback->sol->error ? isl_stat_error : isl_stat_ok;
581 : }
582 :
583 : /* Move down to next level and push callback onto context tableau
584 : * to decrease the level again when it gets rolled back across
585 : * the current state. That is, dec_level will be called with
586 : * the context tableau in the same state as it is when inc_level
587 : * is called.
588 : */
589 0 : static void sol_inc_level(struct isl_sol *sol)
590 : {
591 : struct isl_tab *tab;
592 :
593 0 : if (sol->error)
594 0 : return;
595 :
596 0 : sol->level++;
597 0 : tab = sol->context->op->peek_tab(sol->context);
598 0 : if (isl_tab_push_callback(tab, &sol->dec_level.callback) < 0)
599 0 : sol->error = 1;
600 : }
601 :
602 0 : static void scale_rows(struct isl_mat *mat, isl_int m, int n_row)
603 : {
604 : int i;
605 :
606 0 : if (isl_int_is_one(m))
607 0 : return;
608 :
609 0 : for (i = 0; i < n_row; ++i)
610 0 : isl_seq_scale(mat->row[i], mat->row[i], m, mat->n_col);
611 : }
612 :
613 : /* Add the solution identified by the tableau and the context tableau.
614 : *
615 : * The layout of the variables is as follows.
616 : * tab->n_var is equal to the total number of variables in the input
617 : * map (including divs that were copied from the context)
618 : * + the number of extra divs constructed
619 : * Of these, the first tab->n_param and the last tab->n_div variables
620 : * correspond to the variables in the context, i.e.,
621 : * tab->n_param + tab->n_div = context_tab->n_var
622 : * tab->n_param is equal to the number of parameters and input
623 : * dimensions in the input map
624 : * tab->n_div is equal to the number of divs in the context
625 : *
626 : * If there is no solution, then call add_empty with a basic set
627 : * that corresponds to the context tableau. (If add_empty is NULL,
628 : * then do nothing).
629 : *
630 : * If there is a solution, then first construct a matrix that maps
631 : * all dimensions of the context to the output variables, i.e.,
632 : * the output dimensions in the input map.
633 : * The divs in the input map (if any) that do not correspond to any
634 : * div in the context do not appear in the solution.
635 : * The algorithm will make sure that they have an integer value,
636 : * but these values themselves are of no interest.
637 : * We have to be careful not to drop or rearrange any divs in the
638 : * context because that would change the meaning of the matrix.
639 : *
640 : * To extract the value of the output variables, it should be noted
641 : * that we always use a big parameter M in the main tableau and so
642 : * the variable stored in this tableau is not an output variable x itself, but
643 : * x' = M + x (in case of minimization)
644 : * or
645 : * x' = M - x (in case of maximization)
646 : * If x' appears in a column, then its optimal value is zero,
647 : * which means that the optimal value of x is an unbounded number
648 : * (-M for minimization and M for maximization).
649 : * We currently assume that the output dimensions in the original map
650 : * are bounded, so this cannot occur.
651 : * Similarly, when x' appears in a row, then the coefficient of M in that
652 : * row is necessarily 1.
653 : * If the row in the tableau represents
654 : * d x' = c + d M + e(y)
655 : * then, in case of minimization, the corresponding row in the matrix
656 : * will be
657 : * a c + a e(y)
658 : * with a d = m, the (updated) common denominator of the matrix.
659 : * In case of maximization, the row will be
660 : * -a c - a e(y)
661 : */
662 0 : static void sol_add(struct isl_sol *sol, struct isl_tab *tab)
663 : {
664 0 : struct isl_basic_set *bset = NULL;
665 0 : struct isl_mat *mat = NULL;
666 : unsigned off;
667 : int row;
668 : isl_int m;
669 :
670 0 : if (sol->error || !tab)
671 : goto error;
672 :
673 0 : if (tab->empty && !sol->add_empty)
674 0 : return;
675 0 : if (sol->context->op->is_empty(sol->context))
676 0 : return;
677 :
678 0 : bset = sol_domain(sol);
679 :
680 0 : if (tab->empty) {
681 0 : sol_push_sol(sol, bset, NULL);
682 0 : return;
683 : }
684 :
685 0 : off = 2 + tab->M;
686 :
687 0 : mat = isl_mat_alloc(tab->mat->ctx, 1 + sol->n_out,
688 0 : 1 + tab->n_param + tab->n_div);
689 0 : if (!mat)
690 0 : goto error;
691 :
692 0 : isl_int_init(m);
693 :
694 0 : isl_seq_clr(mat->row[0] + 1, mat->n_col - 1);
695 0 : isl_int_set_si(mat->row[0][0], 1);
696 0 : for (row = 0; row < sol->n_out; ++row) {
697 0 : int i = tab->n_param + row;
698 : int r, j;
699 :
700 0 : isl_seq_clr(mat->row[1 + row], mat->n_col);
701 0 : if (!tab->var[i].is_row) {
702 0 : if (tab->M)
703 0 : isl_die(mat->ctx, isl_error_invalid,
704 : "unbounded optimum", goto error2);
705 0 : continue;
706 : }
707 :
708 0 : r = tab->var[i].index;
709 0 : if (tab->M &&
710 0 : isl_int_ne(tab->mat->row[r][2], tab->mat->row[r][0]))
711 0 : isl_die(mat->ctx, isl_error_invalid,
712 : "unbounded optimum", goto error2);
713 0 : isl_int_gcd(m, mat->row[0][0], tab->mat->row[r][0]);
714 0 : isl_int_divexact(m, tab->mat->row[r][0], m);
715 0 : scale_rows(mat, m, 1 + row);
716 0 : isl_int_divexact(m, mat->row[0][0], tab->mat->row[r][0]);
717 0 : isl_int_mul(mat->row[1 + row][0], m, tab->mat->row[r][1]);
718 0 : for (j = 0; j < tab->n_param; ++j) {
719 : int col;
720 0 : if (tab->var[j].is_row)
721 0 : continue;
722 0 : col = tab->var[j].index;
723 0 : isl_int_mul(mat->row[1 + row][1 + j], m,
724 : tab->mat->row[r][off + col]);
725 : }
726 0 : for (j = 0; j < tab->n_div; ++j) {
727 : int col;
728 0 : if (tab->var[tab->n_var - tab->n_div+j].is_row)
729 0 : continue;
730 0 : col = tab->var[tab->n_var - tab->n_div+j].index;
731 0 : isl_int_mul(mat->row[1 + row][1 + tab->n_param + j], m,
732 : tab->mat->row[r][off + col]);
733 : }
734 0 : if (sol->max)
735 0 : isl_seq_neg(mat->row[1 + row], mat->row[1 + row],
736 : mat->n_col);
737 : }
738 :
739 0 : isl_int_clear(m);
740 :
741 0 : sol_push_sol_mat(sol, bset, mat);
742 0 : return;
743 : error2:
744 0 : isl_int_clear(m);
745 : error:
746 0 : isl_basic_set_free(bset);
747 0 : isl_mat_free(mat);
748 0 : sol->error = 1;
749 : }
750 :
751 : struct isl_sol_map {
752 : struct isl_sol sol;
753 : struct isl_map *map;
754 : struct isl_set *empty;
755 : };
756 :
757 0 : static void sol_map_free(struct isl_sol *sol)
758 : {
759 0 : struct isl_sol_map *sol_map = (struct isl_sol_map *) sol;
760 0 : isl_map_free(sol_map->map);
761 0 : isl_set_free(sol_map->empty);
762 0 : }
763 :
764 : /* This function is called for parts of the context where there is
765 : * no solution, with "bset" corresponding to the context tableau.
766 : * Simply add the basic set to the set "empty".
767 : */
768 0 : static void sol_map_add_empty(struct isl_sol_map *sol,
769 : struct isl_basic_set *bset)
770 : {
771 0 : if (!bset || !sol->empty)
772 : goto error;
773 :
774 0 : sol->empty = isl_set_grow(sol->empty, 1);
775 0 : bset = isl_basic_set_simplify(bset);
776 0 : bset = isl_basic_set_finalize(bset);
777 0 : sol->empty = isl_set_add_basic_set(sol->empty, isl_basic_set_copy(bset));
778 0 : if (!sol->empty)
779 0 : goto error;
780 0 : isl_basic_set_free(bset);
781 0 : return;
782 : error:
783 0 : isl_basic_set_free(bset);
784 0 : sol->sol.error = 1;
785 : }
786 :
787 0 : static void sol_map_add_empty_wrap(struct isl_sol *sol,
788 : struct isl_basic_set *bset)
789 : {
790 0 : sol_map_add_empty((struct isl_sol_map *)sol, bset);
791 0 : }
792 :
793 : /* Given a basic set "dom" that represents the context and a tuple of
794 : * affine expressions "ma" defined over this domain, construct a basic map
795 : * that expresses this function on the domain.
796 : */
797 0 : static void sol_map_add(struct isl_sol_map *sol,
798 : __isl_take isl_basic_set *dom, __isl_take isl_multi_aff *ma)
799 : {
800 : isl_basic_map *bmap;
801 :
802 0 : if (sol->sol.error || !dom || !ma)
803 : goto error;
804 :
805 0 : bmap = isl_basic_map_from_multi_aff2(ma, sol->sol.rational);
806 0 : bmap = isl_basic_map_intersect_domain(bmap, dom);
807 0 : sol->map = isl_map_grow(sol->map, 1);
808 0 : sol->map = isl_map_add_basic_map(sol->map, bmap);
809 0 : if (!sol->map)
810 0 : sol->sol.error = 1;
811 0 : return;
812 : error:
813 0 : isl_basic_set_free(dom);
814 0 : isl_multi_aff_free(ma);
815 0 : sol->sol.error = 1;
816 : }
817 :
818 0 : static void sol_map_add_wrap(struct isl_sol *sol,
819 : __isl_take isl_basic_set *dom, __isl_take isl_multi_aff *ma)
820 : {
821 0 : sol_map_add((struct isl_sol_map *)sol, dom, ma);
822 0 : }
823 :
824 :
825 : /* Store the "parametric constant" of row "row" of tableau "tab" in "line",
826 : * i.e., the constant term and the coefficients of all variables that
827 : * appear in the context tableau.
828 : * Note that the coefficient of the big parameter M is NOT copied.
829 : * The context tableau may not have a big parameter and even when it
830 : * does, it is a different big parameter.
831 : */
832 0 : static void get_row_parameter_line(struct isl_tab *tab, int row, isl_int *line)
833 : {
834 : int i;
835 0 : unsigned off = 2 + tab->M;
836 :
837 0 : isl_int_set(line[0], tab->mat->row[row][1]);
838 0 : for (i = 0; i < tab->n_param; ++i) {
839 0 : if (tab->var[i].is_row)
840 0 : isl_int_set_si(line[1 + i], 0);
841 : else {
842 0 : int col = tab->var[i].index;
843 0 : isl_int_set(line[1 + i], tab->mat->row[row][off + col]);
844 : }
845 : }
846 0 : for (i = 0; i < tab->n_div; ++i) {
847 0 : if (tab->var[tab->n_var - tab->n_div + i].is_row)
848 0 : isl_int_set_si(line[1 + tab->n_param + i], 0);
849 : else {
850 0 : int col = tab->var[tab->n_var - tab->n_div + i].index;
851 0 : isl_int_set(line[1 + tab->n_param + i],
852 : tab->mat->row[row][off + col]);
853 : }
854 : }
855 0 : }
856 :
857 : /* Check if rows "row1" and "row2" have identical "parametric constants",
858 : * as explained above.
859 : * In this case, we also insist that the coefficients of the big parameter
860 : * be the same as the values of the constants will only be the same
861 : * if these coefficients are also the same.
862 : */
863 0 : static int identical_parameter_line(struct isl_tab *tab, int row1, int row2)
864 : {
865 : int i;
866 0 : unsigned off = 2 + tab->M;
867 :
868 0 : if (isl_int_ne(tab->mat->row[row1][1], tab->mat->row[row2][1]))
869 0 : return 0;
870 :
871 0 : if (tab->M && isl_int_ne(tab->mat->row[row1][2],
872 : tab->mat->row[row2][2]))
873 0 : return 0;
874 :
875 0 : for (i = 0; i < tab->n_param + tab->n_div; ++i) {
876 0 : int pos = i < tab->n_param ? i :
877 0 : tab->n_var - tab->n_div + i - tab->n_param;
878 : int col;
879 :
880 0 : if (tab->var[pos].is_row)
881 0 : continue;
882 0 : col = tab->var[pos].index;
883 0 : if (isl_int_ne(tab->mat->row[row1][off + col],
884 : tab->mat->row[row2][off + col]))
885 0 : return 0;
886 : }
887 0 : return 1;
888 : }
889 :
890 : /* Return an inequality that expresses that the "parametric constant"
891 : * should be non-negative.
892 : * This function is only called when the coefficient of the big parameter
893 : * is equal to zero.
894 : */
895 0 : static struct isl_vec *get_row_parameter_ineq(struct isl_tab *tab, int row)
896 : {
897 : struct isl_vec *ineq;
898 :
899 0 : ineq = isl_vec_alloc(tab->mat->ctx, 1 + tab->n_param + tab->n_div);
900 0 : if (!ineq)
901 0 : return NULL;
902 :
903 0 : get_row_parameter_line(tab, row, ineq->el);
904 0 : if (ineq)
905 0 : ineq = isl_vec_normalize(ineq);
906 :
907 0 : return ineq;
908 : }
909 :
910 : /* Normalize a div expression of the form
911 : *
912 : * [(g*f(x) + c)/(g * m)]
913 : *
914 : * with c the constant term and f(x) the remaining coefficients, to
915 : *
916 : * [(f(x) + [c/g])/m]
917 : */
918 0 : static void normalize_div(__isl_keep isl_vec *div)
919 : {
920 0 : isl_ctx *ctx = isl_vec_get_ctx(div);
921 0 : int len = div->size - 2;
922 :
923 0 : isl_seq_gcd(div->el + 2, len, &ctx->normalize_gcd);
924 0 : isl_int_gcd(ctx->normalize_gcd, ctx->normalize_gcd, div->el[0]);
925 :
926 0 : if (isl_int_is_one(ctx->normalize_gcd))
927 0 : return;
928 :
929 0 : isl_int_divexact(div->el[0], div->el[0], ctx->normalize_gcd);
930 0 : isl_int_fdiv_q(div->el[1], div->el[1], ctx->normalize_gcd);
931 0 : isl_seq_scale_down(div->el + 2, div->el + 2, ctx->normalize_gcd, len);
932 : }
933 :
934 : /* Return an integer division for use in a parametric cut based
935 : * on the given row.
936 : * In particular, let the parametric constant of the row be
937 : *
938 : * \sum_i a_i y_i
939 : *
940 : * where y_0 = 1, but none of the y_i corresponds to the big parameter M.
941 : * The div returned is equal to
942 : *
943 : * floor(\sum_i {-a_i} y_i) = floor((\sum_i (-a_i mod d) y_i)/d)
944 : */
945 0 : static struct isl_vec *get_row_parameter_div(struct isl_tab *tab, int row)
946 : {
947 : struct isl_vec *div;
948 :
949 0 : div = isl_vec_alloc(tab->mat->ctx, 1 + 1 + tab->n_param + tab->n_div);
950 0 : if (!div)
951 0 : return NULL;
952 :
953 0 : isl_int_set(div->el[0], tab->mat->row[row][0]);
954 0 : get_row_parameter_line(tab, row, div->el + 1);
955 0 : isl_seq_neg(div->el + 1, div->el + 1, div->size - 1);
956 0 : normalize_div(div);
957 0 : isl_seq_fdiv_r(div->el + 1, div->el + 1, div->el[0], div->size - 1);
958 :
959 0 : return div;
960 : }
961 :
962 : /* Return an integer division for use in transferring an integrality constraint
963 : * to the context.
964 : * In particular, let the parametric constant of the row be
965 : *
966 : * \sum_i a_i y_i
967 : *
968 : * where y_0 = 1, but none of the y_i corresponds to the big parameter M.
969 : * The the returned div is equal to
970 : *
971 : * floor(\sum_i {a_i} y_i) = floor((\sum_i (a_i mod d) y_i)/d)
972 : */
973 0 : static struct isl_vec *get_row_split_div(struct isl_tab *tab, int row)
974 : {
975 : struct isl_vec *div;
976 :
977 0 : div = isl_vec_alloc(tab->mat->ctx, 1 + 1 + tab->n_param + tab->n_div);
978 0 : if (!div)
979 0 : return NULL;
980 :
981 0 : isl_int_set(div->el[0], tab->mat->row[row][0]);
982 0 : get_row_parameter_line(tab, row, div->el + 1);
983 0 : normalize_div(div);
984 0 : isl_seq_fdiv_r(div->el + 1, div->el + 1, div->el[0], div->size - 1);
985 :
986 0 : return div;
987 : }
988 :
989 : /* Construct and return an inequality that expresses an upper bound
990 : * on the given div.
991 : * In particular, if the div is given by
992 : *
993 : * d = floor(e/m)
994 : *
995 : * then the inequality expresses
996 : *
997 : * m d <= e
998 : */
999 0 : static __isl_give isl_vec *ineq_for_div(__isl_keep isl_basic_set *bset,
1000 : unsigned div)
1001 : {
1002 : unsigned total;
1003 : unsigned div_pos;
1004 : struct isl_vec *ineq;
1005 :
1006 0 : if (!bset)
1007 0 : return NULL;
1008 :
1009 0 : total = isl_basic_set_total_dim(bset);
1010 0 : div_pos = 1 + total - bset->n_div + div;
1011 :
1012 0 : ineq = isl_vec_alloc(bset->ctx, 1 + total);
1013 0 : if (!ineq)
1014 0 : return NULL;
1015 :
1016 0 : isl_seq_cpy(ineq->el, bset->div[div] + 1, 1 + total);
1017 0 : isl_int_neg(ineq->el[div_pos], bset->div[div][0]);
1018 0 : return ineq;
1019 : }
1020 :
1021 : /* Given a row in the tableau and a div that was created
1022 : * using get_row_split_div and that has been constrained to equality, i.e.,
1023 : *
1024 : * d = floor(\sum_i {a_i} y_i) = \sum_i {a_i} y_i
1025 : *
1026 : * replace the expression "\sum_i {a_i} y_i" in the row by d,
1027 : * i.e., we subtract "\sum_i {a_i} y_i" and add 1 d.
1028 : * The coefficients of the non-parameters in the tableau have been
1029 : * verified to be integral. We can therefore simply replace coefficient b
1030 : * by floor(b). For the coefficients of the parameters we have
1031 : * floor(a_i) = a_i - {a_i}, while for the other coefficients, we have
1032 : * floor(b) = b.
1033 : */
1034 0 : static struct isl_tab *set_row_cst_to_div(struct isl_tab *tab, int row, int div)
1035 : {
1036 0 : isl_seq_fdiv_q(tab->mat->row[row] + 1, tab->mat->row[row] + 1,
1037 0 : tab->mat->row[row][0], 1 + tab->M + tab->n_col);
1038 :
1039 0 : isl_int_set_si(tab->mat->row[row][0], 1);
1040 :
1041 0 : if (tab->var[tab->n_var - tab->n_div + div].is_row) {
1042 0 : int drow = tab->var[tab->n_var - tab->n_div + div].index;
1043 :
1044 0 : isl_assert(tab->mat->ctx,
1045 : isl_int_is_one(tab->mat->row[drow][0]), goto error);
1046 0 : isl_seq_combine(tab->mat->row[row] + 1,
1047 0 : tab->mat->ctx->one, tab->mat->row[row] + 1,
1048 0 : tab->mat->ctx->one, tab->mat->row[drow] + 1,
1049 0 : 1 + tab->M + tab->n_col);
1050 : } else {
1051 0 : int dcol = tab->var[tab->n_var - tab->n_div + div].index;
1052 :
1053 0 : isl_int_add_ui(tab->mat->row[row][2 + tab->M + dcol],
1054 : tab->mat->row[row][2 + tab->M + dcol], 1);
1055 : }
1056 :
1057 0 : return tab;
1058 : error:
1059 0 : isl_tab_free(tab);
1060 0 : return NULL;
1061 : }
1062 :
1063 : /* Check if the (parametric) constant of the given row is obviously
1064 : * negative, meaning that we don't need to consult the context tableau.
1065 : * If there is a big parameter and its coefficient is non-zero,
1066 : * then this coefficient determines the outcome.
1067 : * Otherwise, we check whether the constant is negative and
1068 : * all non-zero coefficients of parameters are negative and
1069 : * belong to non-negative parameters.
1070 : */
1071 0 : static int is_obviously_neg(struct isl_tab *tab, int row)
1072 : {
1073 : int i;
1074 : int col;
1075 0 : unsigned off = 2 + tab->M;
1076 :
1077 0 : if (tab->M) {
1078 0 : if (isl_int_is_pos(tab->mat->row[row][2]))
1079 0 : return 0;
1080 0 : if (isl_int_is_neg(tab->mat->row[row][2]))
1081 0 : return 1;
1082 : }
1083 :
1084 0 : if (isl_int_is_nonneg(tab->mat->row[row][1]))
1085 0 : return 0;
1086 0 : for (i = 0; i < tab->n_param; ++i) {
1087 : /* Eliminated parameter */
1088 0 : if (tab->var[i].is_row)
1089 0 : continue;
1090 0 : col = tab->var[i].index;
1091 0 : if (isl_int_is_zero(tab->mat->row[row][off + col]))
1092 0 : continue;
1093 0 : if (!tab->var[i].is_nonneg)
1094 0 : return 0;
1095 0 : if (isl_int_is_pos(tab->mat->row[row][off + col]))
1096 0 : return 0;
1097 : }
1098 0 : for (i = 0; i < tab->n_div; ++i) {
1099 0 : if (tab->var[tab->n_var - tab->n_div + i].is_row)
1100 0 : continue;
1101 0 : col = tab->var[tab->n_var - tab->n_div + i].index;
1102 0 : if (isl_int_is_zero(tab->mat->row[row][off + col]))
1103 0 : continue;
1104 0 : if (!tab->var[tab->n_var - tab->n_div + i].is_nonneg)
1105 0 : return 0;
1106 0 : if (isl_int_is_pos(tab->mat->row[row][off + col]))
1107 0 : return 0;
1108 : }
1109 0 : return 1;
1110 : }
1111 :
1112 : /* Check if the (parametric) constant of the given row is obviously
1113 : * non-negative, meaning that we don't need to consult the context tableau.
1114 : * If there is a big parameter and its coefficient is non-zero,
1115 : * then this coefficient determines the outcome.
1116 : * Otherwise, we check whether the constant is non-negative and
1117 : * all non-zero coefficients of parameters are positive and
1118 : * belong to non-negative parameters.
1119 : */
1120 0 : static int is_obviously_nonneg(struct isl_tab *tab, int row)
1121 : {
1122 : int i;
1123 : int col;
1124 0 : unsigned off = 2 + tab->M;
1125 :
1126 0 : if (tab->M) {
1127 0 : if (isl_int_is_pos(tab->mat->row[row][2]))
1128 0 : return 1;
1129 0 : if (isl_int_is_neg(tab->mat->row[row][2]))
1130 0 : return 0;
1131 : }
1132 :
1133 0 : if (isl_int_is_neg(tab->mat->row[row][1]))
1134 0 : return 0;
1135 0 : for (i = 0; i < tab->n_param; ++i) {
1136 : /* Eliminated parameter */
1137 0 : if (tab->var[i].is_row)
1138 0 : continue;
1139 0 : col = tab->var[i].index;
1140 0 : if (isl_int_is_zero(tab->mat->row[row][off + col]))
1141 0 : continue;
1142 0 : if (!tab->var[i].is_nonneg)
1143 0 : return 0;
1144 0 : if (isl_int_is_neg(tab->mat->row[row][off + col]))
1145 0 : return 0;
1146 : }
1147 0 : for (i = 0; i < tab->n_div; ++i) {
1148 0 : if (tab->var[tab->n_var - tab->n_div + i].is_row)
1149 0 : continue;
1150 0 : col = tab->var[tab->n_var - tab->n_div + i].index;
1151 0 : if (isl_int_is_zero(tab->mat->row[row][off + col]))
1152 0 : continue;
1153 0 : if (!tab->var[tab->n_var - tab->n_div + i].is_nonneg)
1154 0 : return 0;
1155 0 : if (isl_int_is_neg(tab->mat->row[row][off + col]))
1156 0 : return 0;
1157 : }
1158 0 : return 1;
1159 : }
1160 :
1161 : /* Given a row r and two columns, return the column that would
1162 : * lead to the lexicographically smallest increment in the sample
1163 : * solution when leaving the basis in favor of the row.
1164 : * Pivoting with column c will increment the sample value by a non-negative
1165 : * constant times a_{V,c}/a_{r,c}, with a_{V,c} the elements of column c
1166 : * corresponding to the non-parametric variables.
1167 : * If variable v appears in a column c_v, then a_{v,c} = 1 iff c = c_v,
1168 : * with all other entries in this virtual row equal to zero.
1169 : * If variable v appears in a row, then a_{v,c} is the element in column c
1170 : * of that row.
1171 : *
1172 : * Let v be the first variable with a_{v,c1}/a_{r,c1} != a_{v,c2}/a_{r,c2}.
1173 : * Then if a_{v,c1}/a_{r,c1} < a_{v,c2}/a_{r,c2}, i.e.,
1174 : * a_{v,c2} a_{r,c1} - a_{v,c1} a_{r,c2} > 0, c1 results in the minimal
1175 : * increment. Otherwise, it's c2.
1176 : */
1177 0 : static int lexmin_col_pair(struct isl_tab *tab,
1178 : int row, int col1, int col2, isl_int tmp)
1179 : {
1180 : int i;
1181 : isl_int *tr;
1182 :
1183 0 : tr = tab->mat->row[row] + 2 + tab->M;
1184 :
1185 0 : for (i = tab->n_param; i < tab->n_var - tab->n_div; ++i) {
1186 : int s1, s2;
1187 : isl_int *r;
1188 :
1189 0 : if (!tab->var[i].is_row) {
1190 0 : if (tab->var[i].index == col1)
1191 0 : return col2;
1192 0 : if (tab->var[i].index == col2)
1193 0 : return col1;
1194 0 : continue;
1195 : }
1196 :
1197 0 : if (tab->var[i].index == row)
1198 0 : continue;
1199 :
1200 0 : r = tab->mat->row[tab->var[i].index] + 2 + tab->M;
1201 0 : s1 = isl_int_sgn(r[col1]);
1202 0 : s2 = isl_int_sgn(r[col2]);
1203 0 : if (s1 == 0 && s2 == 0)
1204 0 : continue;
1205 0 : if (s1 < s2)
1206 0 : return col1;
1207 0 : if (s2 < s1)
1208 0 : return col2;
1209 :
1210 0 : isl_int_mul(tmp, r[col2], tr[col1]);
1211 0 : isl_int_submul(tmp, r[col1], tr[col2]);
1212 0 : if (isl_int_is_pos(tmp))
1213 0 : return col1;
1214 0 : if (isl_int_is_neg(tmp))
1215 0 : return col2;
1216 : }
1217 0 : return -1;
1218 : }
1219 :
1220 : /* Does the index into the tab->var or tab->con array "index"
1221 : * correspond to a variable in the context tableau?
1222 : * In particular, it needs to be an index into the tab->var array and
1223 : * it needs to refer to either one of the first tab->n_param variables or
1224 : * one of the last tab->n_div variables.
1225 : */
1226 0 : static int is_parameter_var(struct isl_tab *tab, int index)
1227 : {
1228 0 : if (index < 0)
1229 0 : return 0;
1230 0 : if (index < tab->n_param)
1231 0 : return 1;
1232 0 : if (index >= tab->n_var - tab->n_div)
1233 0 : return 1;
1234 0 : return 0;
1235 : }
1236 :
1237 : /* Does column "col" of "tab" refer to a variable in the context tableau?
1238 : */
1239 0 : static int col_is_parameter_var(struct isl_tab *tab, int col)
1240 : {
1241 0 : return is_parameter_var(tab, tab->col_var[col]);
1242 : }
1243 :
1244 : /* Does row "row" of "tab" refer to a variable in the context tableau?
1245 : */
1246 0 : static int row_is_parameter_var(struct isl_tab *tab, int row)
1247 : {
1248 0 : return is_parameter_var(tab, tab->row_var[row]);
1249 : }
1250 :
1251 : /* Given a row in the tableau, find and return the column that would
1252 : * result in the lexicographically smallest, but positive, increment
1253 : * in the sample point.
1254 : * If there is no such column, then return tab->n_col.
1255 : * If anything goes wrong, return -1.
1256 : */
1257 0 : static int lexmin_pivot_col(struct isl_tab *tab, int row)
1258 : {
1259 : int j;
1260 0 : int col = tab->n_col;
1261 : isl_int *tr;
1262 : isl_int tmp;
1263 :
1264 0 : tr = tab->mat->row[row] + 2 + tab->M;
1265 :
1266 0 : isl_int_init(tmp);
1267 :
1268 0 : for (j = tab->n_dead; j < tab->n_col; ++j) {
1269 0 : if (col_is_parameter_var(tab, j))
1270 0 : continue;
1271 :
1272 0 : if (!isl_int_is_pos(tr[j]))
1273 0 : continue;
1274 :
1275 0 : if (col == tab->n_col)
1276 0 : col = j;
1277 : else
1278 0 : col = lexmin_col_pair(tab, row, col, j, tmp);
1279 0 : isl_assert(tab->mat->ctx, col >= 0, goto error);
1280 : }
1281 :
1282 0 : isl_int_clear(tmp);
1283 0 : return col;
1284 : error:
1285 0 : isl_int_clear(tmp);
1286 0 : return -1;
1287 : }
1288 :
1289 : /* Return the first known violated constraint, i.e., a non-negative
1290 : * constraint that currently has an either obviously negative value
1291 : * or a previously determined to be negative value.
1292 : *
1293 : * If any constraint has a negative coefficient for the big parameter,
1294 : * if any, then we return one of these first.
1295 : */
1296 0 : static int first_neg(struct isl_tab *tab)
1297 : {
1298 : int row;
1299 :
1300 0 : if (tab->M)
1301 0 : for (row = tab->n_redundant; row < tab->n_row; ++row) {
1302 0 : if (!isl_tab_var_from_row(tab, row)->is_nonneg)
1303 0 : continue;
1304 0 : if (!isl_int_is_neg(tab->mat->row[row][2]))
1305 0 : continue;
1306 0 : if (tab->row_sign)
1307 0 : tab->row_sign[row] = isl_tab_row_neg;
1308 0 : return row;
1309 : }
1310 0 : for (row = tab->n_redundant; row < tab->n_row; ++row) {
1311 0 : if (!isl_tab_var_from_row(tab, row)->is_nonneg)
1312 0 : continue;
1313 0 : if (tab->row_sign) {
1314 0 : if (tab->row_sign[row] == 0 &&
1315 0 : is_obviously_neg(tab, row))
1316 0 : tab->row_sign[row] = isl_tab_row_neg;
1317 0 : if (tab->row_sign[row] != isl_tab_row_neg)
1318 0 : continue;
1319 0 : } else if (!is_obviously_neg(tab, row))
1320 0 : continue;
1321 0 : return row;
1322 : }
1323 0 : return -1;
1324 : }
1325 :
1326 : /* Check whether the invariant that all columns are lexico-positive
1327 : * is satisfied. This function is not called from the current code
1328 : * but is useful during debugging.
1329 : */
1330 : static void check_lexpos(struct isl_tab *tab) __attribute__ ((unused));
1331 0 : static void check_lexpos(struct isl_tab *tab)
1332 : {
1333 0 : unsigned off = 2 + tab->M;
1334 : int col;
1335 : int var;
1336 : int row;
1337 :
1338 0 : for (col = tab->n_dead; col < tab->n_col; ++col) {
1339 0 : if (col_is_parameter_var(tab, col))
1340 0 : continue;
1341 0 : for (var = tab->n_param; var < tab->n_var - tab->n_div; ++var) {
1342 0 : if (!tab->var[var].is_row) {
1343 0 : if (tab->var[var].index == col)
1344 0 : break;
1345 : else
1346 0 : continue;
1347 : }
1348 0 : row = tab->var[var].index;
1349 0 : if (isl_int_is_zero(tab->mat->row[row][off + col]))
1350 0 : continue;
1351 0 : if (isl_int_is_pos(tab->mat->row[row][off + col]))
1352 0 : break;
1353 0 : fprintf(stderr, "lexneg column %d (row %d)\n",
1354 : col, row);
1355 : }
1356 0 : if (var >= tab->n_var - tab->n_div)
1357 0 : fprintf(stderr, "zero column %d\n", col);
1358 : }
1359 0 : }
1360 :
1361 : /* Report to the caller that the given constraint is part of an encountered
1362 : * conflict.
1363 : */
1364 0 : static int report_conflicting_constraint(struct isl_tab *tab, int con)
1365 : {
1366 0 : return tab->conflict(con, tab->conflict_user);
1367 : }
1368 :
1369 : /* Given a conflicting row in the tableau, report all constraints
1370 : * involved in the row to the caller. That is, the row itself
1371 : * (if it represents a constraint) and all constraint columns with
1372 : * non-zero (and therefore negative) coefficients.
1373 : */
1374 0 : static int report_conflict(struct isl_tab *tab, int row)
1375 : {
1376 : int j;
1377 : isl_int *tr;
1378 :
1379 0 : if (!tab->conflict)
1380 0 : return 0;
1381 :
1382 0 : if (tab->row_var[row] < 0 &&
1383 0 : report_conflicting_constraint(tab, ~tab->row_var[row]) < 0)
1384 0 : return -1;
1385 :
1386 0 : tr = tab->mat->row[row] + 2 + tab->M;
1387 :
1388 0 : for (j = tab->n_dead; j < tab->n_col; ++j) {
1389 0 : if (col_is_parameter_var(tab, j))
1390 0 : continue;
1391 :
1392 0 : if (!isl_int_is_neg(tr[j]))
1393 0 : continue;
1394 :
1395 0 : if (tab->col_var[j] < 0 &&
1396 0 : report_conflicting_constraint(tab, ~tab->col_var[j]) < 0)
1397 0 : return -1;
1398 : }
1399 :
1400 0 : return 0;
1401 : }
1402 :
1403 : /* Resolve all known or obviously violated constraints through pivoting.
1404 : * In particular, as long as we can find any violated constraint, we
1405 : * look for a pivoting column that would result in the lexicographically
1406 : * smallest increment in the sample point. If there is no such column
1407 : * then the tableau is infeasible.
1408 : */
1409 : static int restore_lexmin(struct isl_tab *tab) WARN_UNUSED;
1410 0 : static int restore_lexmin(struct isl_tab *tab)
1411 : {
1412 : int row, col;
1413 :
1414 0 : if (!tab)
1415 0 : return -1;
1416 0 : if (tab->empty)
1417 0 : return 0;
1418 0 : while ((row = first_neg(tab)) != -1) {
1419 0 : col = lexmin_pivot_col(tab, row);
1420 0 : if (col >= tab->n_col) {
1421 0 : if (report_conflict(tab, row) < 0)
1422 0 : return -1;
1423 0 : if (isl_tab_mark_empty(tab) < 0)
1424 0 : return -1;
1425 0 : return 0;
1426 : }
1427 0 : if (col < 0)
1428 0 : return -1;
1429 0 : if (isl_tab_pivot(tab, row, col) < 0)
1430 0 : return -1;
1431 : }
1432 0 : return 0;
1433 : }
1434 :
1435 : /* Given a row that represents an equality, look for an appropriate
1436 : * pivoting column.
1437 : * In particular, if there are any non-zero coefficients among
1438 : * the non-parameter variables, then we take the last of these
1439 : * variables. Eliminating this variable in terms of the other
1440 : * variables and/or parameters does not influence the property
1441 : * that all column in the initial tableau are lexicographically
1442 : * positive. The row corresponding to the eliminated variable
1443 : * will only have non-zero entries below the diagonal of the
1444 : * initial tableau. That is, we transform
1445 : *
1446 : * I I
1447 : * 1 into a
1448 : * I I
1449 : *
1450 : * If there is no such non-parameter variable, then we are dealing with
1451 : * pure parameter equality and we pick any parameter with coefficient 1 or -1
1452 : * for elimination. This will ensure that the eliminated parameter
1453 : * always has an integer value whenever all the other parameters are integral.
1454 : * If there is no such parameter then we return -1.
1455 : */
1456 0 : static int last_var_col_or_int_par_col(struct isl_tab *tab, int row)
1457 : {
1458 0 : unsigned off = 2 + tab->M;
1459 : int i;
1460 :
1461 0 : for (i = tab->n_var - tab->n_div - 1; i >= 0 && i >= tab->n_param; --i) {
1462 : int col;
1463 0 : if (tab->var[i].is_row)
1464 0 : continue;
1465 0 : col = tab->var[i].index;
1466 0 : if (col <= tab->n_dead)
1467 0 : continue;
1468 0 : if (!isl_int_is_zero(tab->mat->row[row][off + col]))
1469 0 : return col;
1470 : }
1471 0 : for (i = tab->n_dead; i < tab->n_col; ++i) {
1472 0 : if (isl_int_is_one(tab->mat->row[row][off + i]))
1473 0 : return i;
1474 0 : if (isl_int_is_negone(tab->mat->row[row][off + i]))
1475 0 : return i;
1476 : }
1477 0 : return -1;
1478 : }
1479 :
1480 : /* Add an equality that is known to be valid to the tableau.
1481 : * We first check if we can eliminate a variable or a parameter.
1482 : * If not, we add the equality as two inequalities.
1483 : * In this case, the equality was a pure parameter equality and there
1484 : * is no need to resolve any constraint violations.
1485 : *
1486 : * This function assumes that at least two more rows and at least
1487 : * two more elements in the constraint array are available in the tableau.
1488 : */
1489 0 : static struct isl_tab *add_lexmin_valid_eq(struct isl_tab *tab, isl_int *eq)
1490 : {
1491 : int i;
1492 : int r;
1493 :
1494 0 : if (!tab)
1495 0 : return NULL;
1496 0 : r = isl_tab_add_row(tab, eq);
1497 0 : if (r < 0)
1498 0 : goto error;
1499 :
1500 0 : r = tab->con[r].index;
1501 0 : i = last_var_col_or_int_par_col(tab, r);
1502 0 : if (i < 0) {
1503 0 : tab->con[r].is_nonneg = 1;
1504 0 : if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r]) < 0)
1505 0 : goto error;
1506 0 : isl_seq_neg(eq, eq, 1 + tab->n_var);
1507 0 : r = isl_tab_add_row(tab, eq);
1508 0 : if (r < 0)
1509 0 : goto error;
1510 0 : tab->con[r].is_nonneg = 1;
1511 0 : if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r]) < 0)
1512 0 : goto error;
1513 : } else {
1514 0 : if (isl_tab_pivot(tab, r, i) < 0)
1515 0 : goto error;
1516 0 : if (isl_tab_kill_col(tab, i) < 0)
1517 0 : goto error;
1518 0 : tab->n_eq++;
1519 : }
1520 :
1521 0 : return tab;
1522 : error:
1523 0 : isl_tab_free(tab);
1524 0 : return NULL;
1525 : }
1526 :
1527 : /* Check if the given row is a pure constant.
1528 : */
1529 0 : static int is_constant(struct isl_tab *tab, int row)
1530 : {
1531 0 : unsigned off = 2 + tab->M;
1532 :
1533 0 : return isl_seq_first_non_zero(tab->mat->row[row] + off + tab->n_dead,
1534 0 : tab->n_col - tab->n_dead) == -1;
1535 : }
1536 :
1537 : /* Is the given row a parametric constant?
1538 : * That is, does it only involve variables that also appear in the context?
1539 : */
1540 0 : static int is_parametric_constant(struct isl_tab *tab, int row)
1541 : {
1542 0 : unsigned off = 2 + tab->M;
1543 : int col;
1544 :
1545 0 : for (col = tab->n_dead; col < tab->n_col; ++col) {
1546 0 : if (col_is_parameter_var(tab, col))
1547 0 : continue;
1548 0 : if (isl_int_is_zero(tab->mat->row[row][off + col]))
1549 0 : continue;
1550 0 : return 0;
1551 : }
1552 :
1553 0 : return 1;
1554 : }
1555 :
1556 : /* Add an equality that may or may not be valid to the tableau.
1557 : * If the resulting row is a pure constant, then it must be zero.
1558 : * Otherwise, the resulting tableau is empty.
1559 : *
1560 : * If the row is not a pure constant, then we add two inequalities,
1561 : * each time checking that they can be satisfied.
1562 : * In the end we try to use one of the two constraints to eliminate
1563 : * a column.
1564 : *
1565 : * This function assumes that at least two more rows and at least
1566 : * two more elements in the constraint array are available in the tableau.
1567 : */
1568 : static int add_lexmin_eq(struct isl_tab *tab, isl_int *eq) WARN_UNUSED;
1569 0 : static int add_lexmin_eq(struct isl_tab *tab, isl_int *eq)
1570 : {
1571 : int r1, r2;
1572 : int row;
1573 : struct isl_tab_undo *snap;
1574 :
1575 0 : if (!tab)
1576 0 : return -1;
1577 0 : snap = isl_tab_snap(tab);
1578 0 : r1 = isl_tab_add_row(tab, eq);
1579 0 : if (r1 < 0)
1580 0 : return -1;
1581 0 : tab->con[r1].is_nonneg = 1;
1582 0 : if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r1]) < 0)
1583 0 : return -1;
1584 :
1585 0 : row = tab->con[r1].index;
1586 0 : if (is_constant(tab, row)) {
1587 0 : if (!isl_int_is_zero(tab->mat->row[row][1]) ||
1588 0 : (tab->M && !isl_int_is_zero(tab->mat->row[row][2]))) {
1589 0 : if (isl_tab_mark_empty(tab) < 0)
1590 0 : return -1;
1591 0 : return 0;
1592 : }
1593 0 : if (isl_tab_rollback(tab, snap) < 0)
1594 0 : return -1;
1595 0 : return 0;
1596 : }
1597 :
1598 0 : if (restore_lexmin(tab) < 0)
1599 0 : return -1;
1600 0 : if (tab->empty)
1601 0 : return 0;
1602 :
1603 0 : isl_seq_neg(eq, eq, 1 + tab->n_var);
1604 :
1605 0 : r2 = isl_tab_add_row(tab, eq);
1606 0 : if (r2 < 0)
1607 0 : return -1;
1608 0 : tab->con[r2].is_nonneg = 1;
1609 0 : if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r2]) < 0)
1610 0 : return -1;
1611 :
1612 0 : if (restore_lexmin(tab) < 0)
1613 0 : return -1;
1614 0 : if (tab->empty)
1615 0 : return 0;
1616 :
1617 0 : if (!tab->con[r1].is_row) {
1618 0 : if (isl_tab_kill_col(tab, tab->con[r1].index) < 0)
1619 0 : return -1;
1620 0 : } else if (!tab->con[r2].is_row) {
1621 0 : if (isl_tab_kill_col(tab, tab->con[r2].index) < 0)
1622 0 : return -1;
1623 : }
1624 :
1625 0 : if (tab->bmap) {
1626 0 : tab->bmap = isl_basic_map_add_ineq(tab->bmap, eq);
1627 0 : if (isl_tab_push(tab, isl_tab_undo_bmap_ineq) < 0)
1628 0 : return -1;
1629 0 : isl_seq_neg(eq, eq, 1 + tab->n_var);
1630 0 : tab->bmap = isl_basic_map_add_ineq(tab->bmap, eq);
1631 0 : isl_seq_neg(eq, eq, 1 + tab->n_var);
1632 0 : if (isl_tab_push(tab, isl_tab_undo_bmap_ineq) < 0)
1633 0 : return -1;
1634 0 : if (!tab->bmap)
1635 0 : return -1;
1636 : }
1637 :
1638 0 : return 0;
1639 : }
1640 :
1641 : /* Add an inequality to the tableau, resolving violations using
1642 : * restore_lexmin.
1643 : *
1644 : * This function assumes that at least one more row and at least
1645 : * one more element in the constraint array are available in the tableau.
1646 : */
1647 0 : static struct isl_tab *add_lexmin_ineq(struct isl_tab *tab, isl_int *ineq)
1648 : {
1649 : int r;
1650 :
1651 0 : if (!tab)
1652 0 : return NULL;
1653 0 : if (tab->bmap) {
1654 0 : tab->bmap = isl_basic_map_add_ineq(tab->bmap, ineq);
1655 0 : if (isl_tab_push(tab, isl_tab_undo_bmap_ineq) < 0)
1656 0 : goto error;
1657 0 : if (!tab->bmap)
1658 0 : goto error;
1659 : }
1660 0 : r = isl_tab_add_row(tab, ineq);
1661 0 : if (r < 0)
1662 0 : goto error;
1663 0 : tab->con[r].is_nonneg = 1;
1664 0 : if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r]) < 0)
1665 0 : goto error;
1666 0 : if (isl_tab_row_is_redundant(tab, tab->con[r].index)) {
1667 0 : if (isl_tab_mark_redundant(tab, tab->con[r].index) < 0)
1668 0 : goto error;
1669 0 : return tab;
1670 : }
1671 :
1672 0 : if (restore_lexmin(tab) < 0)
1673 0 : goto error;
1674 0 : if (!tab->empty && tab->con[r].is_row &&
1675 0 : isl_tab_row_is_redundant(tab, tab->con[r].index))
1676 0 : if (isl_tab_mark_redundant(tab, tab->con[r].index) < 0)
1677 0 : goto error;
1678 0 : return tab;
1679 : error:
1680 0 : isl_tab_free(tab);
1681 0 : return NULL;
1682 : }
1683 :
1684 : /* Check if the coefficients of the parameters are all integral.
1685 : */
1686 0 : static int integer_parameter(struct isl_tab *tab, int row)
1687 : {
1688 : int i;
1689 : int col;
1690 0 : unsigned off = 2 + tab->M;
1691 :
1692 0 : for (i = 0; i < tab->n_param; ++i) {
1693 : /* Eliminated parameter */
1694 0 : if (tab->var[i].is_row)
1695 0 : continue;
1696 0 : col = tab->var[i].index;
1697 0 : if (!isl_int_is_divisible_by(tab->mat->row[row][off + col],
1698 : tab->mat->row[row][0]))
1699 0 : return 0;
1700 : }
1701 0 : for (i = 0; i < tab->n_div; ++i) {
1702 0 : if (tab->var[tab->n_var - tab->n_div + i].is_row)
1703 0 : continue;
1704 0 : col = tab->var[tab->n_var - tab->n_div + i].index;
1705 0 : if (!isl_int_is_divisible_by(tab->mat->row[row][off + col],
1706 : tab->mat->row[row][0]))
1707 0 : return 0;
1708 : }
1709 0 : return 1;
1710 : }
1711 :
1712 : /* Check if the coefficients of the non-parameter variables are all integral.
1713 : */
1714 0 : static int integer_variable(struct isl_tab *tab, int row)
1715 : {
1716 : int i;
1717 0 : unsigned off = 2 + tab->M;
1718 :
1719 0 : for (i = tab->n_dead; i < tab->n_col; ++i) {
1720 0 : if (col_is_parameter_var(tab, i))
1721 0 : continue;
1722 0 : if (!isl_int_is_divisible_by(tab->mat->row[row][off + i],
1723 : tab->mat->row[row][0]))
1724 0 : return 0;
1725 : }
1726 0 : return 1;
1727 : }
1728 :
1729 : /* Check if the constant term is integral.
1730 : */
1731 0 : static int integer_constant(struct isl_tab *tab, int row)
1732 : {
1733 0 : return isl_int_is_divisible_by(tab->mat->row[row][1],
1734 : tab->mat->row[row][0]);
1735 : }
1736 :
1737 : #define I_CST 1 << 0
1738 : #define I_PAR 1 << 1
1739 : #define I_VAR 1 << 2
1740 :
1741 : /* Check for next (non-parameter) variable after "var" (first if var == -1)
1742 : * that is non-integer and therefore requires a cut and return
1743 : * the index of the variable.
1744 : * For parametric tableaus, there are three parts in a row,
1745 : * the constant, the coefficients of the parameters and the rest.
1746 : * For each part, we check whether the coefficients in that part
1747 : * are all integral and if so, set the corresponding flag in *f.
1748 : * If the constant and the parameter part are integral, then the
1749 : * current sample value is integral and no cut is required
1750 : * (irrespective of whether the variable part is integral).
1751 : */
1752 0 : static int next_non_integer_var(struct isl_tab *tab, int var, int *f)
1753 : {
1754 0 : var = var < 0 ? tab->n_param : var + 1;
1755 :
1756 0 : for (; var < tab->n_var - tab->n_div; ++var) {
1757 0 : int flags = 0;
1758 : int row;
1759 0 : if (!tab->var[var].is_row)
1760 0 : continue;
1761 0 : row = tab->var[var].index;
1762 0 : if (integer_constant(tab, row))
1763 0 : ISL_FL_SET(flags, I_CST);
1764 0 : if (integer_parameter(tab, row))
1765 0 : ISL_FL_SET(flags, I_PAR);
1766 0 : if (ISL_FL_ISSET(flags, I_CST) && ISL_FL_ISSET(flags, I_PAR))
1767 0 : continue;
1768 0 : if (integer_variable(tab, row))
1769 0 : ISL_FL_SET(flags, I_VAR);
1770 0 : *f = flags;
1771 0 : return var;
1772 : }
1773 0 : return -1;
1774 : }
1775 :
1776 : /* Check for first (non-parameter) variable that is non-integer and
1777 : * therefore requires a cut and return the corresponding row.
1778 : * For parametric tableaus, there are three parts in a row,
1779 : * the constant, the coefficients of the parameters and the rest.
1780 : * For each part, we check whether the coefficients in that part
1781 : * are all integral and if so, set the corresponding flag in *f.
1782 : * If the constant and the parameter part are integral, then the
1783 : * current sample value is integral and no cut is required
1784 : * (irrespective of whether the variable part is integral).
1785 : */
1786 0 : static int first_non_integer_row(struct isl_tab *tab, int *f)
1787 : {
1788 0 : int var = next_non_integer_var(tab, -1, f);
1789 :
1790 0 : return var < 0 ? -1 : tab->var[var].index;
1791 : }
1792 :
1793 : /* Add a (non-parametric) cut to cut away the non-integral sample
1794 : * value of the given row.
1795 : *
1796 : * If the row is given by
1797 : *
1798 : * m r = f + \sum_i a_i y_i
1799 : *
1800 : * then the cut is
1801 : *
1802 : * c = - {-f/m} + \sum_i {a_i/m} y_i >= 0
1803 : *
1804 : * The big parameter, if any, is ignored, since it is assumed to be big
1805 : * enough to be divisible by any integer.
1806 : * If the tableau is actually a parametric tableau, then this function
1807 : * is only called when all coefficients of the parameters are integral.
1808 : * The cut therefore has zero coefficients for the parameters.
1809 : *
1810 : * The current value is known to be negative, so row_sign, if it
1811 : * exists, is set accordingly.
1812 : *
1813 : * Return the row of the cut or -1.
1814 : */
1815 0 : static int add_cut(struct isl_tab *tab, int row)
1816 : {
1817 : int i;
1818 : int r;
1819 : isl_int *r_row;
1820 0 : unsigned off = 2 + tab->M;
1821 :
1822 0 : if (isl_tab_extend_cons(tab, 1) < 0)
1823 0 : return -1;
1824 0 : r = isl_tab_allocate_con(tab);
1825 0 : if (r < 0)
1826 0 : return -1;
1827 :
1828 0 : r_row = tab->mat->row[tab->con[r].index];
1829 0 : isl_int_set(r_row[0], tab->mat->row[row][0]);
1830 0 : isl_int_neg(r_row[1], tab->mat->row[row][1]);
1831 0 : isl_int_fdiv_r(r_row[1], r_row[1], tab->mat->row[row][0]);
1832 0 : isl_int_neg(r_row[1], r_row[1]);
1833 0 : if (tab->M)
1834 0 : isl_int_set_si(r_row[2], 0);
1835 0 : for (i = 0; i < tab->n_col; ++i)
1836 0 : isl_int_fdiv_r(r_row[off + i],
1837 : tab->mat->row[row][off + i], tab->mat->row[row][0]);
1838 :
1839 0 : tab->con[r].is_nonneg = 1;
1840 0 : if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r]) < 0)
1841 0 : return -1;
1842 0 : if (tab->row_sign)
1843 0 : tab->row_sign[tab->con[r].index] = isl_tab_row_neg;
1844 :
1845 0 : return tab->con[r].index;
1846 : }
1847 :
1848 : #define CUT_ALL 1
1849 : #define CUT_ONE 0
1850 :
1851 : /* Given a non-parametric tableau, add cuts until an integer
1852 : * sample point is obtained or until the tableau is determined
1853 : * to be integer infeasible.
1854 : * As long as there is any non-integer value in the sample point,
1855 : * we add appropriate cuts, if possible, for each of these
1856 : * non-integer values and then resolve the violated
1857 : * cut constraints using restore_lexmin.
1858 : * If one of the corresponding rows is equal to an integral
1859 : * combination of variables/constraints plus a non-integral constant,
1860 : * then there is no way to obtain an integer point and we return
1861 : * a tableau that is marked empty.
1862 : * The parameter cutting_strategy controls the strategy used when adding cuts
1863 : * to remove non-integer points. CUT_ALL adds all possible cuts
1864 : * before continuing the search. CUT_ONE adds only one cut at a time.
1865 : */
1866 0 : static struct isl_tab *cut_to_integer_lexmin(struct isl_tab *tab,
1867 : int cutting_strategy)
1868 : {
1869 : int var;
1870 : int row;
1871 : int flags;
1872 :
1873 0 : if (!tab)
1874 0 : return NULL;
1875 0 : if (tab->empty)
1876 0 : return tab;
1877 :
1878 0 : while ((var = next_non_integer_var(tab, -1, &flags)) != -1) {
1879 : do {
1880 0 : if (ISL_FL_ISSET(flags, I_VAR)) {
1881 0 : if (isl_tab_mark_empty(tab) < 0)
1882 0 : goto error;
1883 0 : return tab;
1884 : }
1885 0 : row = tab->var[var].index;
1886 0 : row = add_cut(tab, row);
1887 0 : if (row < 0)
1888 0 : goto error;
1889 0 : if (cutting_strategy == CUT_ONE)
1890 0 : break;
1891 0 : } while ((var = next_non_integer_var(tab, var, &flags)) != -1);
1892 0 : if (restore_lexmin(tab) < 0)
1893 0 : goto error;
1894 0 : if (tab->empty)
1895 0 : break;
1896 : }
1897 0 : return tab;
1898 : error:
1899 0 : isl_tab_free(tab);
1900 0 : return NULL;
1901 : }
1902 :
1903 : /* Check whether all the currently active samples also satisfy the inequality
1904 : * "ineq" (treated as an equality if eq is set).
1905 : * Remove those samples that do not.
1906 : */
1907 0 : static struct isl_tab *check_samples(struct isl_tab *tab, isl_int *ineq, int eq)
1908 : {
1909 : int i;
1910 : isl_int v;
1911 :
1912 0 : if (!tab)
1913 0 : return NULL;
1914 :
1915 0 : isl_assert(tab->mat->ctx, tab->bmap, goto error);
1916 0 : isl_assert(tab->mat->ctx, tab->samples, goto error);
1917 0 : isl_assert(tab->mat->ctx, tab->samples->n_col == 1 + tab->n_var, goto error);
1918 :
1919 0 : isl_int_init(v);
1920 0 : for (i = tab->n_outside; i < tab->n_sample; ++i) {
1921 : int sgn;
1922 0 : isl_seq_inner_product(ineq, tab->samples->row[i],
1923 0 : 1 + tab->n_var, &v);
1924 0 : sgn = isl_int_sgn(v);
1925 0 : if (eq ? (sgn == 0) : (sgn >= 0))
1926 0 : continue;
1927 0 : tab = isl_tab_drop_sample(tab, i);
1928 0 : if (!tab)
1929 0 : break;
1930 : }
1931 0 : isl_int_clear(v);
1932 :
1933 0 : return tab;
1934 : error:
1935 0 : isl_tab_free(tab);
1936 0 : return NULL;
1937 : }
1938 :
1939 : /* Check whether the sample value of the tableau is finite,
1940 : * i.e., either the tableau does not use a big parameter, or
1941 : * all values of the variables are equal to the big parameter plus
1942 : * some constant. This constant is the actual sample value.
1943 : */
1944 0 : static int sample_is_finite(struct isl_tab *tab)
1945 : {
1946 : int i;
1947 :
1948 0 : if (!tab->M)
1949 0 : return 1;
1950 :
1951 0 : for (i = 0; i < tab->n_var; ++i) {
1952 : int row;
1953 0 : if (!tab->var[i].is_row)
1954 0 : return 0;
1955 0 : row = tab->var[i].index;
1956 0 : if (isl_int_ne(tab->mat->row[row][0], tab->mat->row[row][2]))
1957 0 : return 0;
1958 : }
1959 0 : return 1;
1960 : }
1961 :
1962 : /* Check if the context tableau of sol has any integer points.
1963 : * Leave tab in empty state if no integer point can be found.
1964 : * If an integer point can be found and if moreover it is finite,
1965 : * then it is added to the list of sample values.
1966 : *
1967 : * This function is only called when none of the currently active sample
1968 : * values satisfies the most recently added constraint.
1969 : */
1970 0 : static struct isl_tab *check_integer_feasible(struct isl_tab *tab)
1971 : {
1972 : struct isl_tab_undo *snap;
1973 :
1974 0 : if (!tab)
1975 0 : return NULL;
1976 :
1977 0 : snap = isl_tab_snap(tab);
1978 0 : if (isl_tab_push_basis(tab) < 0)
1979 0 : goto error;
1980 :
1981 0 : tab = cut_to_integer_lexmin(tab, CUT_ALL);
1982 0 : if (!tab)
1983 0 : goto error;
1984 :
1985 0 : if (!tab->empty && sample_is_finite(tab)) {
1986 : struct isl_vec *sample;
1987 :
1988 0 : sample = isl_tab_get_sample_value(tab);
1989 :
1990 0 : if (isl_tab_add_sample(tab, sample) < 0)
1991 0 : goto error;
1992 : }
1993 :
1994 0 : if (!tab->empty && isl_tab_rollback(tab, snap) < 0)
1995 0 : goto error;
1996 :
1997 0 : return tab;
1998 : error:
1999 0 : isl_tab_free(tab);
2000 0 : return NULL;
2001 : }
2002 :
2003 : /* Check if any of the currently active sample values satisfies
2004 : * the inequality "ineq" (an equality if eq is set).
2005 : */
2006 0 : static int tab_has_valid_sample(struct isl_tab *tab, isl_int *ineq, int eq)
2007 : {
2008 : int i;
2009 : isl_int v;
2010 :
2011 0 : if (!tab)
2012 0 : return -1;
2013 :
2014 0 : isl_assert(tab->mat->ctx, tab->bmap, return -1);
2015 0 : isl_assert(tab->mat->ctx, tab->samples, return -1);
2016 0 : isl_assert(tab->mat->ctx, tab->samples->n_col == 1 + tab->n_var, return -1);
2017 :
2018 0 : isl_int_init(v);
2019 0 : for (i = tab->n_outside; i < tab->n_sample; ++i) {
2020 : int sgn;
2021 0 : isl_seq_inner_product(ineq, tab->samples->row[i],
2022 0 : 1 + tab->n_var, &v);
2023 0 : sgn = isl_int_sgn(v);
2024 0 : if (eq ? (sgn == 0) : (sgn >= 0))
2025 0 : break;
2026 : }
2027 0 : isl_int_clear(v);
2028 :
2029 0 : return i < tab->n_sample;
2030 : }
2031 :
2032 : /* Insert a div specified by "div" to the tableau "tab" at position "pos" and
2033 : * return isl_bool_true if the div is obviously non-negative.
2034 : */
2035 0 : static isl_bool context_tab_insert_div(struct isl_tab *tab, int pos,
2036 : __isl_keep isl_vec *div,
2037 : isl_stat (*add_ineq)(void *user, isl_int *), void *user)
2038 : {
2039 : int i;
2040 : int r;
2041 : struct isl_mat *samples;
2042 : int nonneg;
2043 :
2044 0 : r = isl_tab_insert_div(tab, pos, div, add_ineq, user);
2045 0 : if (r < 0)
2046 0 : return isl_bool_error;
2047 0 : nonneg = tab->var[r].is_nonneg;
2048 0 : tab->var[r].frozen = 1;
2049 :
2050 0 : samples = isl_mat_extend(tab->samples,
2051 0 : tab->n_sample, 1 + tab->n_var);
2052 0 : tab->samples = samples;
2053 0 : if (!samples)
2054 0 : return isl_bool_error;
2055 0 : for (i = tab->n_outside; i < samples->n_row; ++i) {
2056 0 : isl_seq_inner_product(div->el + 1, samples->row[i],
2057 0 : div->size - 1, &samples->row[i][samples->n_col - 1]);
2058 0 : isl_int_fdiv_q(samples->row[i][samples->n_col - 1],
2059 : samples->row[i][samples->n_col - 1], div->el[0]);
2060 : }
2061 0 : tab->samples = isl_mat_move_cols(tab->samples, 1 + pos,
2062 : 1 + tab->n_var - 1, 1);
2063 0 : if (!tab->samples)
2064 0 : return isl_bool_error;
2065 :
2066 0 : return nonneg;
2067 : }
2068 :
2069 : /* Add a div specified by "div" to both the main tableau and
2070 : * the context tableau. In case of the main tableau, we only
2071 : * need to add an extra div. In the context tableau, we also
2072 : * need to express the meaning of the div.
2073 : * Return the index of the div or -1 if anything went wrong.
2074 : *
2075 : * The new integer division is added before any unknown integer
2076 : * divisions in the context to ensure that it does not get
2077 : * equated to some linear combination involving unknown integer
2078 : * divisions.
2079 : */
2080 0 : static int add_div(struct isl_tab *tab, struct isl_context *context,
2081 : __isl_keep isl_vec *div)
2082 : {
2083 : int r;
2084 : int pos;
2085 : isl_bool nonneg;
2086 0 : struct isl_tab *context_tab = context->op->peek_tab(context);
2087 :
2088 0 : if (!tab || !context_tab)
2089 : goto error;
2090 :
2091 0 : pos = context_tab->n_var - context->n_unknown;
2092 0 : if ((nonneg = context->op->insert_div(context, pos, div)) < 0)
2093 0 : goto error;
2094 :
2095 0 : if (!context->op->is_ok(context))
2096 0 : goto error;
2097 :
2098 0 : pos = tab->n_var - context->n_unknown;
2099 0 : if (isl_tab_extend_vars(tab, 1) < 0)
2100 0 : goto error;
2101 0 : r = isl_tab_insert_var(tab, pos);
2102 0 : if (r < 0)
2103 0 : goto error;
2104 0 : if (nonneg)
2105 0 : tab->var[r].is_nonneg = 1;
2106 0 : tab->var[r].frozen = 1;
2107 0 : tab->n_div++;
2108 :
2109 0 : return tab->n_div - 1 - context->n_unknown;
2110 : error:
2111 0 : context->op->invalidate(context);
2112 0 : return -1;
2113 : }
2114 :
2115 0 : static int find_div(struct isl_tab *tab, isl_int *div, isl_int denom)
2116 : {
2117 : int i;
2118 0 : unsigned total = isl_basic_map_total_dim(tab->bmap);
2119 :
2120 0 : for (i = 0; i < tab->bmap->n_div; ++i) {
2121 0 : if (isl_int_ne(tab->bmap->div[i][0], denom))
2122 0 : continue;
2123 0 : if (!isl_seq_eq(tab->bmap->div[i] + 1, div, 1 + total))
2124 0 : continue;
2125 0 : return i;
2126 : }
2127 0 : return -1;
2128 : }
2129 :
2130 : /* Return the index of a div that corresponds to "div".
2131 : * We first check if we already have such a div and if not, we create one.
2132 : */
2133 0 : static int get_div(struct isl_tab *tab, struct isl_context *context,
2134 : struct isl_vec *div)
2135 : {
2136 : int d;
2137 0 : struct isl_tab *context_tab = context->op->peek_tab(context);
2138 :
2139 0 : if (!context_tab)
2140 0 : return -1;
2141 :
2142 0 : d = find_div(context_tab, div->el + 1, div->el[0]);
2143 0 : if (d != -1)
2144 0 : return d;
2145 :
2146 0 : return add_div(tab, context, div);
2147 : }
2148 :
2149 : /* Add a parametric cut to cut away the non-integral sample value
2150 : * of the given row.
2151 : * Let a_i be the coefficients of the constant term and the parameters
2152 : * and let b_i be the coefficients of the variables or constraints
2153 : * in basis of the tableau.
2154 : * Let q be the div q = floor(\sum_i {-a_i} y_i).
2155 : *
2156 : * The cut is expressed as
2157 : *
2158 : * c = \sum_i -{-a_i} y_i + \sum_i {b_i} x_i + q >= 0
2159 : *
2160 : * If q did not already exist in the context tableau, then it is added first.
2161 : * If q is in a column of the main tableau then the "+ q" can be accomplished
2162 : * by setting the corresponding entry to the denominator of the constraint.
2163 : * If q happens to be in a row of the main tableau, then the corresponding
2164 : * row needs to be added instead (taking care of the denominators).
2165 : * Note that this is very unlikely, but perhaps not entirely impossible.
2166 : *
2167 : * The current value of the cut is known to be negative (or at least
2168 : * non-positive), so row_sign is set accordingly.
2169 : *
2170 : * Return the row of the cut or -1.
2171 : */
2172 0 : static int add_parametric_cut(struct isl_tab *tab, int row,
2173 : struct isl_context *context)
2174 : {
2175 : struct isl_vec *div;
2176 : int d;
2177 : int i;
2178 : int r;
2179 : isl_int *r_row;
2180 : int col;
2181 : int n;
2182 0 : unsigned off = 2 + tab->M;
2183 :
2184 0 : if (!context)
2185 0 : return -1;
2186 :
2187 0 : div = get_row_parameter_div(tab, row);
2188 0 : if (!div)
2189 0 : return -1;
2190 :
2191 0 : n = tab->n_div - context->n_unknown;
2192 0 : d = context->op->get_div(context, tab, div);
2193 0 : isl_vec_free(div);
2194 0 : if (d < 0)
2195 0 : return -1;
2196 :
2197 0 : if (isl_tab_extend_cons(tab, 1) < 0)
2198 0 : return -1;
2199 0 : r = isl_tab_allocate_con(tab);
2200 0 : if (r < 0)
2201 0 : return -1;
2202 :
2203 0 : r_row = tab->mat->row[tab->con[r].index];
2204 0 : isl_int_set(r_row[0], tab->mat->row[row][0]);
2205 0 : isl_int_neg(r_row[1], tab->mat->row[row][1]);
2206 0 : isl_int_fdiv_r(r_row[1], r_row[1], tab->mat->row[row][0]);
2207 0 : isl_int_neg(r_row[1], r_row[1]);
2208 0 : if (tab->M)
2209 0 : isl_int_set_si(r_row[2], 0);
2210 0 : for (i = 0; i < tab->n_param; ++i) {
2211 0 : if (tab->var[i].is_row)
2212 0 : continue;
2213 0 : col = tab->var[i].index;
2214 0 : isl_int_neg(r_row[off + col], tab->mat->row[row][off + col]);
2215 0 : isl_int_fdiv_r(r_row[off + col], r_row[off + col],
2216 : tab->mat->row[row][0]);
2217 0 : isl_int_neg(r_row[off + col], r_row[off + col]);
2218 : }
2219 0 : for (i = 0; i < tab->n_div; ++i) {
2220 0 : if (tab->var[tab->n_var - tab->n_div + i].is_row)
2221 0 : continue;
2222 0 : col = tab->var[tab->n_var - tab->n_div + i].index;
2223 0 : isl_int_neg(r_row[off + col], tab->mat->row[row][off + col]);
2224 0 : isl_int_fdiv_r(r_row[off + col], r_row[off + col],
2225 : tab->mat->row[row][0]);
2226 0 : isl_int_neg(r_row[off + col], r_row[off + col]);
2227 : }
2228 0 : for (i = 0; i < tab->n_col; ++i) {
2229 0 : if (tab->col_var[i] >= 0 &&
2230 0 : (tab->col_var[i] < tab->n_param ||
2231 0 : tab->col_var[i] >= tab->n_var - tab->n_div))
2232 0 : continue;
2233 0 : isl_int_fdiv_r(r_row[off + i],
2234 : tab->mat->row[row][off + i], tab->mat->row[row][0]);
2235 : }
2236 0 : if (tab->var[tab->n_var - tab->n_div + d].is_row) {
2237 : isl_int gcd;
2238 0 : int d_row = tab->var[tab->n_var - tab->n_div + d].index;
2239 0 : isl_int_init(gcd);
2240 0 : isl_int_gcd(gcd, tab->mat->row[d_row][0], r_row[0]);
2241 0 : isl_int_divexact(r_row[0], r_row[0], gcd);
2242 0 : isl_int_divexact(gcd, tab->mat->row[d_row][0], gcd);
2243 0 : isl_seq_combine(r_row + 1, gcd, r_row + 1,
2244 0 : r_row[0], tab->mat->row[d_row] + 1,
2245 0 : off - 1 + tab->n_col);
2246 0 : isl_int_mul(r_row[0], r_row[0], tab->mat->row[d_row][0]);
2247 0 : isl_int_clear(gcd);
2248 : } else {
2249 0 : col = tab->var[tab->n_var - tab->n_div + d].index;
2250 0 : isl_int_set(r_row[off + col], tab->mat->row[row][0]);
2251 : }
2252 :
2253 0 : tab->con[r].is_nonneg = 1;
2254 0 : if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r]) < 0)
2255 0 : return -1;
2256 0 : if (tab->row_sign)
2257 0 : tab->row_sign[tab->con[r].index] = isl_tab_row_neg;
2258 :
2259 0 : row = tab->con[r].index;
2260 :
2261 0 : if (d >= n && context->op->detect_equalities(context, tab) < 0)
2262 0 : return -1;
2263 :
2264 0 : return row;
2265 : }
2266 :
2267 : /* Construct a tableau for bmap that can be used for computing
2268 : * the lexicographic minimum (or maximum) of bmap.
2269 : * If not NULL, then dom is the domain where the minimum
2270 : * should be computed. In this case, we set up a parametric
2271 : * tableau with row signs (initialized to "unknown").
2272 : * If M is set, then the tableau will use a big parameter.
2273 : * If max is set, then a maximum should be computed instead of a minimum.
2274 : * This means that for each variable x, the tableau will contain the variable
2275 : * x' = M - x, rather than x' = M + x. This in turn means that the coefficient
2276 : * of the variables in all constraints are negated prior to adding them
2277 : * to the tableau.
2278 : */
2279 0 : static __isl_give struct isl_tab *tab_for_lexmin(__isl_keep isl_basic_map *bmap,
2280 : __isl_keep isl_basic_set *dom, unsigned M, int max)
2281 : {
2282 : int i;
2283 : struct isl_tab *tab;
2284 : unsigned n_var;
2285 : unsigned o_var;
2286 :
2287 0 : tab = isl_tab_alloc(bmap->ctx, 2 * bmap->n_eq + bmap->n_ineq + 1,
2288 : isl_basic_map_total_dim(bmap), M);
2289 0 : if (!tab)
2290 0 : return NULL;
2291 :
2292 0 : tab->rational = ISL_F_ISSET(bmap, ISL_BASIC_MAP_RATIONAL);
2293 0 : if (dom) {
2294 0 : tab->n_param = isl_basic_set_total_dim(dom) - dom->n_div;
2295 0 : tab->n_div = dom->n_div;
2296 0 : tab->row_sign = isl_calloc_array(bmap->ctx,
2297 : enum isl_tab_row_sign, tab->mat->n_row);
2298 0 : if (tab->mat->n_row && !tab->row_sign)
2299 0 : goto error;
2300 : }
2301 0 : if (ISL_F_ISSET(bmap, ISL_BASIC_MAP_EMPTY)) {
2302 0 : if (isl_tab_mark_empty(tab) < 0)
2303 0 : goto error;
2304 0 : return tab;
2305 : }
2306 :
2307 0 : for (i = tab->n_param; i < tab->n_var - tab->n_div; ++i) {
2308 0 : tab->var[i].is_nonneg = 1;
2309 0 : tab->var[i].frozen = 1;
2310 : }
2311 0 : o_var = 1 + tab->n_param;
2312 0 : n_var = tab->n_var - tab->n_param - tab->n_div;
2313 0 : for (i = 0; i < bmap->n_eq; ++i) {
2314 0 : if (max)
2315 0 : isl_seq_neg(bmap->eq[i] + o_var,
2316 0 : bmap->eq[i] + o_var, n_var);
2317 0 : tab = add_lexmin_valid_eq(tab, bmap->eq[i]);
2318 0 : if (max)
2319 0 : isl_seq_neg(bmap->eq[i] + o_var,
2320 0 : bmap->eq[i] + o_var, n_var);
2321 0 : if (!tab || tab->empty)
2322 0 : return tab;
2323 : }
2324 0 : if (bmap->n_eq && restore_lexmin(tab) < 0)
2325 0 : goto error;
2326 0 : for (i = 0; i < bmap->n_ineq; ++i) {
2327 0 : if (max)
2328 0 : isl_seq_neg(bmap->ineq[i] + o_var,
2329 0 : bmap->ineq[i] + o_var, n_var);
2330 0 : tab = add_lexmin_ineq(tab, bmap->ineq[i]);
2331 0 : if (max)
2332 0 : isl_seq_neg(bmap->ineq[i] + o_var,
2333 0 : bmap->ineq[i] + o_var, n_var);
2334 0 : if (!tab || tab->empty)
2335 0 : return tab;
2336 : }
2337 0 : return tab;
2338 : error:
2339 0 : isl_tab_free(tab);
2340 0 : return NULL;
2341 : }
2342 :
2343 : /* Given a main tableau where more than one row requires a split,
2344 : * determine and return the "best" row to split on.
2345 : *
2346 : * If any of the rows requiring a split only involves
2347 : * variables that also appear in the context tableau,
2348 : * then the negative part is guaranteed not to have a solution.
2349 : * It is therefore best to split on any of these rows first.
2350 : *
2351 : * Otherwise,
2352 : * given two rows in the main tableau, if the inequality corresponding
2353 : * to the first row is redundant with respect to that of the second row
2354 : * in the current tableau, then it is better to split on the second row,
2355 : * since in the positive part, both rows will be positive.
2356 : * (In the negative part a pivot will have to be performed and just about
2357 : * anything can happen to the sign of the other row.)
2358 : *
2359 : * As a simple heuristic, we therefore select the row that makes the most
2360 : * of the other rows redundant.
2361 : *
2362 : * Perhaps it would also be useful to look at the number of constraints
2363 : * that conflict with any given constraint.
2364 : *
2365 : * best is the best row so far (-1 when we have not found any row yet).
2366 : * best_r is the number of other rows made redundant by row best.
2367 : * When best is still -1, bset_r is meaningless, but it is initialized
2368 : * to some arbitrary value (0) anyway. Without this redundant initialization
2369 : * valgrind may warn about uninitialized memory accesses when isl
2370 : * is compiled with some versions of gcc.
2371 : */
2372 0 : static int best_split(struct isl_tab *tab, struct isl_tab *context_tab)
2373 : {
2374 : struct isl_tab_undo *snap;
2375 : int split;
2376 : int row;
2377 0 : int best = -1;
2378 0 : int best_r = 0;
2379 :
2380 0 : if (isl_tab_extend_cons(context_tab, 2) < 0)
2381 0 : return -1;
2382 :
2383 0 : snap = isl_tab_snap(context_tab);
2384 :
2385 0 : for (split = tab->n_redundant; split < tab->n_row; ++split) {
2386 : struct isl_tab_undo *snap2;
2387 0 : struct isl_vec *ineq = NULL;
2388 0 : int r = 0;
2389 : int ok;
2390 :
2391 0 : if (!isl_tab_var_from_row(tab, split)->is_nonneg)
2392 0 : continue;
2393 0 : if (tab->row_sign[split] != isl_tab_row_any)
2394 0 : continue;
2395 :
2396 0 : if (is_parametric_constant(tab, split))
2397 0 : return split;
2398 :
2399 0 : ineq = get_row_parameter_ineq(tab, split);
2400 0 : if (!ineq)
2401 0 : return -1;
2402 0 : ok = isl_tab_add_ineq(context_tab, ineq->el) >= 0;
2403 0 : isl_vec_free(ineq);
2404 0 : if (!ok)
2405 0 : return -1;
2406 :
2407 0 : snap2 = isl_tab_snap(context_tab);
2408 :
2409 0 : for (row = tab->n_redundant; row < tab->n_row; ++row) {
2410 : struct isl_tab_var *var;
2411 :
2412 0 : if (row == split)
2413 0 : continue;
2414 0 : if (!isl_tab_var_from_row(tab, row)->is_nonneg)
2415 0 : continue;
2416 0 : if (tab->row_sign[row] != isl_tab_row_any)
2417 0 : continue;
2418 :
2419 0 : ineq = get_row_parameter_ineq(tab, row);
2420 0 : if (!ineq)
2421 0 : return -1;
2422 0 : ok = isl_tab_add_ineq(context_tab, ineq->el) >= 0;
2423 0 : isl_vec_free(ineq);
2424 0 : if (!ok)
2425 0 : return -1;
2426 0 : var = &context_tab->con[context_tab->n_con - 1];
2427 0 : if (!context_tab->empty &&
2428 0 : !isl_tab_min_at_most_neg_one(context_tab, var))
2429 0 : r++;
2430 0 : if (isl_tab_rollback(context_tab, snap2) < 0)
2431 0 : return -1;
2432 : }
2433 0 : if (best == -1 || r > best_r) {
2434 0 : best = split;
2435 0 : best_r = r;
2436 : }
2437 0 : if (isl_tab_rollback(context_tab, snap) < 0)
2438 0 : return -1;
2439 : }
2440 :
2441 0 : return best;
2442 : }
2443 :
2444 0 : static struct isl_basic_set *context_lex_peek_basic_set(
2445 : struct isl_context *context)
2446 : {
2447 0 : struct isl_context_lex *clex = (struct isl_context_lex *)context;
2448 0 : if (!clex->tab)
2449 0 : return NULL;
2450 0 : return isl_tab_peek_bset(clex->tab);
2451 : }
2452 :
2453 0 : static struct isl_tab *context_lex_peek_tab(struct isl_context *context)
2454 : {
2455 0 : struct isl_context_lex *clex = (struct isl_context_lex *)context;
2456 0 : return clex->tab;
2457 : }
2458 :
2459 0 : static void context_lex_add_eq(struct isl_context *context, isl_int *eq,
2460 : int check, int update)
2461 : {
2462 0 : struct isl_context_lex *clex = (struct isl_context_lex *)context;
2463 0 : if (isl_tab_extend_cons(clex->tab, 2) < 0)
2464 0 : goto error;
2465 0 : if (add_lexmin_eq(clex->tab, eq) < 0)
2466 0 : goto error;
2467 0 : if (check) {
2468 0 : int v = tab_has_valid_sample(clex->tab, eq, 1);
2469 0 : if (v < 0)
2470 0 : goto error;
2471 0 : if (!v)
2472 0 : clex->tab = check_integer_feasible(clex->tab);
2473 : }
2474 0 : if (update)
2475 0 : clex->tab = check_samples(clex->tab, eq, 1);
2476 0 : return;
2477 : error:
2478 0 : isl_tab_free(clex->tab);
2479 0 : clex->tab = NULL;
2480 : }
2481 :
2482 0 : static void context_lex_add_ineq(struct isl_context *context, isl_int *ineq,
2483 : int check, int update)
2484 : {
2485 0 : struct isl_context_lex *clex = (struct isl_context_lex *)context;
2486 0 : if (isl_tab_extend_cons(clex->tab, 1) < 0)
2487 0 : goto error;
2488 0 : clex->tab = add_lexmin_ineq(clex->tab, ineq);
2489 0 : if (check) {
2490 0 : int v = tab_has_valid_sample(clex->tab, ineq, 0);
2491 0 : if (v < 0)
2492 0 : goto error;
2493 0 : if (!v)
2494 0 : clex->tab = check_integer_feasible(clex->tab);
2495 : }
2496 0 : if (update)
2497 0 : clex->tab = check_samples(clex->tab, ineq, 0);
2498 0 : return;
2499 : error:
2500 0 : isl_tab_free(clex->tab);
2501 0 : clex->tab = NULL;
2502 : }
2503 :
2504 0 : static isl_stat context_lex_add_ineq_wrap(void *user, isl_int *ineq)
2505 : {
2506 0 : struct isl_context *context = (struct isl_context *)user;
2507 0 : context_lex_add_ineq(context, ineq, 0, 0);
2508 0 : return context->op->is_ok(context) ? isl_stat_ok : isl_stat_error;
2509 : }
2510 :
2511 : /* Check which signs can be obtained by "ineq" on all the currently
2512 : * active sample values. See row_sign for more information.
2513 : */
2514 0 : static enum isl_tab_row_sign tab_ineq_sign(struct isl_tab *tab, isl_int *ineq,
2515 : int strict)
2516 : {
2517 : int i;
2518 : int sgn;
2519 : isl_int tmp;
2520 0 : enum isl_tab_row_sign res = isl_tab_row_unknown;
2521 :
2522 0 : isl_assert(tab->mat->ctx, tab->samples, return isl_tab_row_unknown);
2523 0 : isl_assert(tab->mat->ctx, tab->samples->n_col == 1 + tab->n_var,
2524 : return isl_tab_row_unknown);
2525 :
2526 0 : isl_int_init(tmp);
2527 0 : for (i = tab->n_outside; i < tab->n_sample; ++i) {
2528 0 : isl_seq_inner_product(tab->samples->row[i], ineq,
2529 0 : 1 + tab->n_var, &tmp);
2530 0 : sgn = isl_int_sgn(tmp);
2531 0 : if (sgn > 0 || (sgn == 0 && strict)) {
2532 0 : if (res == isl_tab_row_unknown)
2533 0 : res = isl_tab_row_pos;
2534 0 : if (res == isl_tab_row_neg)
2535 0 : res = isl_tab_row_any;
2536 : }
2537 0 : if (sgn < 0) {
2538 0 : if (res == isl_tab_row_unknown)
2539 0 : res = isl_tab_row_neg;
2540 0 : if (res == isl_tab_row_pos)
2541 0 : res = isl_tab_row_any;
2542 : }
2543 0 : if (res == isl_tab_row_any)
2544 0 : break;
2545 : }
2546 0 : isl_int_clear(tmp);
2547 :
2548 0 : return res;
2549 : }
2550 :
2551 0 : static enum isl_tab_row_sign context_lex_ineq_sign(struct isl_context *context,
2552 : isl_int *ineq, int strict)
2553 : {
2554 0 : struct isl_context_lex *clex = (struct isl_context_lex *)context;
2555 0 : return tab_ineq_sign(clex->tab, ineq, strict);
2556 : }
2557 :
2558 : /* Check whether "ineq" can be added to the tableau without rendering
2559 : * it infeasible.
2560 : */
2561 0 : static int context_lex_test_ineq(struct isl_context *context, isl_int *ineq)
2562 : {
2563 0 : struct isl_context_lex *clex = (struct isl_context_lex *)context;
2564 : struct isl_tab_undo *snap;
2565 : int feasible;
2566 :
2567 0 : if (!clex->tab)
2568 0 : return -1;
2569 :
2570 0 : if (isl_tab_extend_cons(clex->tab, 1) < 0)
2571 0 : return -1;
2572 :
2573 0 : snap = isl_tab_snap(clex->tab);
2574 0 : if (isl_tab_push_basis(clex->tab) < 0)
2575 0 : return -1;
2576 0 : clex->tab = add_lexmin_ineq(clex->tab, ineq);
2577 0 : clex->tab = check_integer_feasible(clex->tab);
2578 0 : if (!clex->tab)
2579 0 : return -1;
2580 0 : feasible = !clex->tab->empty;
2581 0 : if (isl_tab_rollback(clex->tab, snap) < 0)
2582 0 : return -1;
2583 :
2584 0 : return feasible;
2585 : }
2586 :
2587 0 : static int context_lex_get_div(struct isl_context *context, struct isl_tab *tab,
2588 : struct isl_vec *div)
2589 : {
2590 0 : return get_div(tab, context, div);
2591 : }
2592 :
2593 : /* Insert a div specified by "div" to the context tableau at position "pos" and
2594 : * return isl_bool_true if the div is obviously non-negative.
2595 : * context_tab_add_div will always return isl_bool_true, because all variables
2596 : * in a isl_context_lex tableau are non-negative.
2597 : * However, if we are using a big parameter in the context, then this only
2598 : * reflects the non-negativity of the variable used to _encode_ the
2599 : * div, i.e., div' = M + div, so we can't draw any conclusions.
2600 : */
2601 0 : static isl_bool context_lex_insert_div(struct isl_context *context, int pos,
2602 : __isl_keep isl_vec *div)
2603 : {
2604 0 : struct isl_context_lex *clex = (struct isl_context_lex *)context;
2605 : isl_bool nonneg;
2606 0 : nonneg = context_tab_insert_div(clex->tab, pos, div,
2607 : context_lex_add_ineq_wrap, context);
2608 0 : if (nonneg < 0)
2609 0 : return isl_bool_error;
2610 0 : if (clex->tab->M)
2611 0 : return isl_bool_false;
2612 0 : return nonneg;
2613 : }
2614 :
2615 0 : static int context_lex_detect_equalities(struct isl_context *context,
2616 : struct isl_tab *tab)
2617 : {
2618 0 : return 0;
2619 : }
2620 :
2621 0 : static int context_lex_best_split(struct isl_context *context,
2622 : struct isl_tab *tab)
2623 : {
2624 0 : struct isl_context_lex *clex = (struct isl_context_lex *)context;
2625 : struct isl_tab_undo *snap;
2626 : int r;
2627 :
2628 0 : snap = isl_tab_snap(clex->tab);
2629 0 : if (isl_tab_push_basis(clex->tab) < 0)
2630 0 : return -1;
2631 0 : r = best_split(tab, clex->tab);
2632 :
2633 0 : if (r >= 0 && isl_tab_rollback(clex->tab, snap) < 0)
2634 0 : return -1;
2635 :
2636 0 : return r;
2637 : }
2638 :
2639 0 : static int context_lex_is_empty(struct isl_context *context)
2640 : {
2641 0 : struct isl_context_lex *clex = (struct isl_context_lex *)context;
2642 0 : if (!clex->tab)
2643 0 : return -1;
2644 0 : return clex->tab->empty;
2645 : }
2646 :
2647 0 : static void *context_lex_save(struct isl_context *context)
2648 : {
2649 0 : struct isl_context_lex *clex = (struct isl_context_lex *)context;
2650 : struct isl_tab_undo *snap;
2651 :
2652 0 : snap = isl_tab_snap(clex->tab);
2653 0 : if (isl_tab_push_basis(clex->tab) < 0)
2654 0 : return NULL;
2655 0 : if (isl_tab_save_samples(clex->tab) < 0)
2656 0 : return NULL;
2657 :
2658 0 : return snap;
2659 : }
2660 :
2661 0 : static void context_lex_restore(struct isl_context *context, void *save)
2662 : {
2663 0 : struct isl_context_lex *clex = (struct isl_context_lex *)context;
2664 0 : if (isl_tab_rollback(clex->tab, (struct isl_tab_undo *)save) < 0) {
2665 0 : isl_tab_free(clex->tab);
2666 0 : clex->tab = NULL;
2667 : }
2668 0 : }
2669 :
2670 0 : static void context_lex_discard(void *save)
2671 : {
2672 0 : }
2673 :
2674 0 : static int context_lex_is_ok(struct isl_context *context)
2675 : {
2676 0 : struct isl_context_lex *clex = (struct isl_context_lex *)context;
2677 0 : return !!clex->tab;
2678 : }
2679 :
2680 : /* For each variable in the context tableau, check if the variable can
2681 : * only attain non-negative values. If so, mark the parameter as non-negative
2682 : * in the main tableau. This allows for a more direct identification of some
2683 : * cases of violated constraints.
2684 : */
2685 0 : static struct isl_tab *tab_detect_nonnegative_parameters(struct isl_tab *tab,
2686 : struct isl_tab *context_tab)
2687 : {
2688 : int i;
2689 : struct isl_tab_undo *snap;
2690 0 : struct isl_vec *ineq = NULL;
2691 : struct isl_tab_var *var;
2692 : int n;
2693 :
2694 0 : if (context_tab->n_var == 0)
2695 0 : return tab;
2696 :
2697 0 : ineq = isl_vec_alloc(tab->mat->ctx, 1 + context_tab->n_var);
2698 0 : if (!ineq)
2699 0 : goto error;
2700 :
2701 0 : if (isl_tab_extend_cons(context_tab, 1) < 0)
2702 0 : goto error;
2703 :
2704 0 : snap = isl_tab_snap(context_tab);
2705 :
2706 0 : n = 0;
2707 0 : isl_seq_clr(ineq->el, ineq->size);
2708 0 : for (i = 0; i < context_tab->n_var; ++i) {
2709 0 : isl_int_set_si(ineq->el[1 + i], 1);
2710 0 : if (isl_tab_add_ineq(context_tab, ineq->el) < 0)
2711 0 : goto error;
2712 0 : var = &context_tab->con[context_tab->n_con - 1];
2713 0 : if (!context_tab->empty &&
2714 0 : !isl_tab_min_at_most_neg_one(context_tab, var)) {
2715 0 : int j = i;
2716 0 : if (i >= tab->n_param)
2717 0 : j = i - tab->n_param + tab->n_var - tab->n_div;
2718 0 : tab->var[j].is_nonneg = 1;
2719 0 : n++;
2720 : }
2721 0 : isl_int_set_si(ineq->el[1 + i], 0);
2722 0 : if (isl_tab_rollback(context_tab, snap) < 0)
2723 0 : goto error;
2724 : }
2725 :
2726 0 : if (context_tab->M && n == context_tab->n_var) {
2727 0 : context_tab->mat = isl_mat_drop_cols(context_tab->mat, 2, 1);
2728 0 : context_tab->M = 0;
2729 : }
2730 :
2731 0 : isl_vec_free(ineq);
2732 0 : return tab;
2733 : error:
2734 0 : isl_vec_free(ineq);
2735 0 : isl_tab_free(tab);
2736 0 : return NULL;
2737 : }
2738 :
2739 0 : static struct isl_tab *context_lex_detect_nonnegative_parameters(
2740 : struct isl_context *context, struct isl_tab *tab)
2741 : {
2742 0 : struct isl_context_lex *clex = (struct isl_context_lex *)context;
2743 : struct isl_tab_undo *snap;
2744 :
2745 0 : if (!tab)
2746 0 : return NULL;
2747 :
2748 0 : snap = isl_tab_snap(clex->tab);
2749 0 : if (isl_tab_push_basis(clex->tab) < 0)
2750 0 : goto error;
2751 :
2752 0 : tab = tab_detect_nonnegative_parameters(tab, clex->tab);
2753 :
2754 0 : if (isl_tab_rollback(clex->tab, snap) < 0)
2755 0 : goto error;
2756 :
2757 0 : return tab;
2758 : error:
2759 0 : isl_tab_free(tab);
2760 0 : return NULL;
2761 : }
2762 :
2763 0 : static void context_lex_invalidate(struct isl_context *context)
2764 : {
2765 0 : struct isl_context_lex *clex = (struct isl_context_lex *)context;
2766 0 : isl_tab_free(clex->tab);
2767 0 : clex->tab = NULL;
2768 0 : }
2769 :
2770 0 : static __isl_null struct isl_context *context_lex_free(
2771 : struct isl_context *context)
2772 : {
2773 0 : struct isl_context_lex *clex = (struct isl_context_lex *)context;
2774 0 : isl_tab_free(clex->tab);
2775 0 : free(clex);
2776 :
2777 0 : return NULL;
2778 : }
2779 :
2780 : struct isl_context_op isl_context_lex_op = {
2781 : context_lex_detect_nonnegative_parameters,
2782 : context_lex_peek_basic_set,
2783 : context_lex_peek_tab,
2784 : context_lex_add_eq,
2785 : context_lex_add_ineq,
2786 : context_lex_ineq_sign,
2787 : context_lex_test_ineq,
2788 : context_lex_get_div,
2789 : context_lex_insert_div,
2790 : context_lex_detect_equalities,
2791 : context_lex_best_split,
2792 : context_lex_is_empty,
2793 : context_lex_is_ok,
2794 : context_lex_save,
2795 : context_lex_restore,
2796 : context_lex_discard,
2797 : context_lex_invalidate,
2798 : context_lex_free,
2799 : };
2800 :
2801 0 : static struct isl_tab *context_tab_for_lexmin(__isl_take isl_basic_set *bset)
2802 : {
2803 : struct isl_tab *tab;
2804 :
2805 0 : if (!bset)
2806 0 : return NULL;
2807 0 : tab = tab_for_lexmin(bset_to_bmap(bset), NULL, 1, 0);
2808 0 : if (isl_tab_track_bset(tab, bset) < 0)
2809 0 : goto error;
2810 0 : tab = isl_tab_init_samples(tab);
2811 0 : return tab;
2812 : error:
2813 0 : isl_tab_free(tab);
2814 0 : return NULL;
2815 : }
2816 :
2817 0 : static struct isl_context *isl_context_lex_alloc(struct isl_basic_set *dom)
2818 : {
2819 : struct isl_context_lex *clex;
2820 :
2821 0 : if (!dom)
2822 0 : return NULL;
2823 :
2824 0 : clex = isl_alloc_type(dom->ctx, struct isl_context_lex);
2825 0 : if (!clex)
2826 0 : return NULL;
2827 :
2828 0 : clex->context.op = &isl_context_lex_op;
2829 :
2830 0 : clex->tab = context_tab_for_lexmin(isl_basic_set_copy(dom));
2831 0 : if (restore_lexmin(clex->tab) < 0)
2832 0 : goto error;
2833 0 : clex->tab = check_integer_feasible(clex->tab);
2834 0 : if (!clex->tab)
2835 0 : goto error;
2836 :
2837 0 : return &clex->context;
2838 : error:
2839 0 : clex->context.op->free(&clex->context);
2840 0 : return NULL;
2841 : }
2842 :
2843 : /* Representation of the context when using generalized basis reduction.
2844 : *
2845 : * "shifted" contains the offsets of the unit hypercubes that lie inside the
2846 : * context. Any rational point in "shifted" can therefore be rounded
2847 : * up to an integer point in the context.
2848 : * If the context is constrained by any equality, then "shifted" is not used
2849 : * as it would be empty.
2850 : */
2851 : struct isl_context_gbr {
2852 : struct isl_context context;
2853 : struct isl_tab *tab;
2854 : struct isl_tab *shifted;
2855 : struct isl_tab *cone;
2856 : };
2857 :
2858 0 : static struct isl_tab *context_gbr_detect_nonnegative_parameters(
2859 : struct isl_context *context, struct isl_tab *tab)
2860 : {
2861 0 : struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
2862 0 : if (!tab)
2863 0 : return NULL;
2864 0 : return tab_detect_nonnegative_parameters(tab, cgbr->tab);
2865 : }
2866 :
2867 0 : static struct isl_basic_set *context_gbr_peek_basic_set(
2868 : struct isl_context *context)
2869 : {
2870 0 : struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
2871 0 : if (!cgbr->tab)
2872 0 : return NULL;
2873 0 : return isl_tab_peek_bset(cgbr->tab);
2874 : }
2875 :
2876 0 : static struct isl_tab *context_gbr_peek_tab(struct isl_context *context)
2877 : {
2878 0 : struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
2879 0 : return cgbr->tab;
2880 : }
2881 :
2882 : /* Initialize the "shifted" tableau of the context, which
2883 : * contains the constraints of the original tableau shifted
2884 : * by the sum of all negative coefficients. This ensures
2885 : * that any rational point in the shifted tableau can
2886 : * be rounded up to yield an integer point in the original tableau.
2887 : */
2888 0 : static void gbr_init_shifted(struct isl_context_gbr *cgbr)
2889 : {
2890 : int i, j;
2891 : struct isl_vec *cst;
2892 0 : struct isl_basic_set *bset = isl_tab_peek_bset(cgbr->tab);
2893 0 : unsigned dim = isl_basic_set_total_dim(bset);
2894 :
2895 0 : cst = isl_vec_alloc(cgbr->tab->mat->ctx, bset->n_ineq);
2896 0 : if (!cst)
2897 0 : return;
2898 :
2899 0 : for (i = 0; i < bset->n_ineq; ++i) {
2900 0 : isl_int_set(cst->el[i], bset->ineq[i][0]);
2901 0 : for (j = 0; j < dim; ++j) {
2902 0 : if (!isl_int_is_neg(bset->ineq[i][1 + j]))
2903 0 : continue;
2904 0 : isl_int_add(bset->ineq[i][0], bset->ineq[i][0],
2905 : bset->ineq[i][1 + j]);
2906 : }
2907 : }
2908 :
2909 0 : cgbr->shifted = isl_tab_from_basic_set(bset, 0);
2910 :
2911 0 : for (i = 0; i < bset->n_ineq; ++i)
2912 0 : isl_int_set(bset->ineq[i][0], cst->el[i]);
2913 :
2914 0 : isl_vec_free(cst);
2915 : }
2916 :
2917 : /* Check if the shifted tableau is non-empty, and if so
2918 : * use the sample point to construct an integer point
2919 : * of the context tableau.
2920 : */
2921 0 : static struct isl_vec *gbr_get_shifted_sample(struct isl_context_gbr *cgbr)
2922 : {
2923 : struct isl_vec *sample;
2924 :
2925 0 : if (!cgbr->shifted)
2926 0 : gbr_init_shifted(cgbr);
2927 0 : if (!cgbr->shifted)
2928 0 : return NULL;
2929 0 : if (cgbr->shifted->empty)
2930 0 : return isl_vec_alloc(cgbr->tab->mat->ctx, 0);
2931 :
2932 0 : sample = isl_tab_get_sample_value(cgbr->shifted);
2933 0 : sample = isl_vec_ceil(sample);
2934 :
2935 0 : return sample;
2936 : }
2937 :
2938 0 : static __isl_give isl_basic_set *drop_constant_terms(
2939 : __isl_take isl_basic_set *bset)
2940 : {
2941 : int i;
2942 :
2943 0 : if (!bset)
2944 0 : return NULL;
2945 :
2946 0 : for (i = 0; i < bset->n_eq; ++i)
2947 0 : isl_int_set_si(bset->eq[i][0], 0);
2948 :
2949 0 : for (i = 0; i < bset->n_ineq; ++i)
2950 0 : isl_int_set_si(bset->ineq[i][0], 0);
2951 :
2952 0 : return bset;
2953 : }
2954 :
2955 0 : static int use_shifted(struct isl_context_gbr *cgbr)
2956 : {
2957 0 : if (!cgbr->tab)
2958 0 : return 0;
2959 0 : return cgbr->tab->bmap->n_eq == 0 && cgbr->tab->bmap->n_div == 0;
2960 : }
2961 :
2962 0 : static struct isl_vec *gbr_get_sample(struct isl_context_gbr *cgbr)
2963 : {
2964 : struct isl_basic_set *bset;
2965 : struct isl_basic_set *cone;
2966 :
2967 0 : if (isl_tab_sample_is_integer(cgbr->tab))
2968 0 : return isl_tab_get_sample_value(cgbr->tab);
2969 :
2970 0 : if (use_shifted(cgbr)) {
2971 : struct isl_vec *sample;
2972 :
2973 0 : sample = gbr_get_shifted_sample(cgbr);
2974 0 : if (!sample || sample->size > 0)
2975 0 : return sample;
2976 :
2977 0 : isl_vec_free(sample);
2978 : }
2979 :
2980 0 : if (!cgbr->cone) {
2981 0 : bset = isl_tab_peek_bset(cgbr->tab);
2982 0 : cgbr->cone = isl_tab_from_recession_cone(bset, 0);
2983 0 : if (!cgbr->cone)
2984 0 : return NULL;
2985 0 : if (isl_tab_track_bset(cgbr->cone,
2986 : isl_basic_set_copy(bset)) < 0)
2987 0 : return NULL;
2988 : }
2989 0 : if (isl_tab_detect_implicit_equalities(cgbr->cone) < 0)
2990 0 : return NULL;
2991 :
2992 0 : if (cgbr->cone->n_dead == cgbr->cone->n_col) {
2993 : struct isl_vec *sample;
2994 : struct isl_tab_undo *snap;
2995 :
2996 0 : if (cgbr->tab->basis) {
2997 0 : if (cgbr->tab->basis->n_col != 1 + cgbr->tab->n_var) {
2998 0 : isl_mat_free(cgbr->tab->basis);
2999 0 : cgbr->tab->basis = NULL;
3000 : }
3001 0 : cgbr->tab->n_zero = 0;
3002 0 : cgbr->tab->n_unbounded = 0;
3003 : }
3004 :
3005 0 : snap = isl_tab_snap(cgbr->tab);
3006 :
3007 0 : sample = isl_tab_sample(cgbr->tab);
3008 :
3009 0 : if (!sample || isl_tab_rollback(cgbr->tab, snap) < 0) {
3010 0 : isl_vec_free(sample);
3011 0 : return NULL;
3012 : }
3013 :
3014 0 : return sample;
3015 : }
3016 :
3017 0 : cone = isl_basic_set_dup(isl_tab_peek_bset(cgbr->cone));
3018 0 : cone = drop_constant_terms(cone);
3019 0 : cone = isl_basic_set_update_from_tab(cone, cgbr->cone);
3020 0 : cone = isl_basic_set_underlying_set(cone);
3021 0 : cone = isl_basic_set_gauss(cone, NULL);
3022 :
3023 0 : bset = isl_basic_set_dup(isl_tab_peek_bset(cgbr->tab));
3024 0 : bset = isl_basic_set_update_from_tab(bset, cgbr->tab);
3025 0 : bset = isl_basic_set_underlying_set(bset);
3026 0 : bset = isl_basic_set_gauss(bset, NULL);
3027 :
3028 0 : return isl_basic_set_sample_with_cone(bset, cone);
3029 : }
3030 :
3031 0 : static void check_gbr_integer_feasible(struct isl_context_gbr *cgbr)
3032 : {
3033 : struct isl_vec *sample;
3034 :
3035 0 : if (!cgbr->tab)
3036 0 : return;
3037 :
3038 0 : if (cgbr->tab->empty)
3039 0 : return;
3040 :
3041 0 : sample = gbr_get_sample(cgbr);
3042 0 : if (!sample)
3043 0 : goto error;
3044 :
3045 0 : if (sample->size == 0) {
3046 0 : isl_vec_free(sample);
3047 0 : if (isl_tab_mark_empty(cgbr->tab) < 0)
3048 0 : goto error;
3049 0 : return;
3050 : }
3051 :
3052 0 : if (isl_tab_add_sample(cgbr->tab, sample) < 0)
3053 0 : goto error;
3054 :
3055 0 : return;
3056 : error:
3057 0 : isl_tab_free(cgbr->tab);
3058 0 : cgbr->tab = NULL;
3059 : }
3060 :
3061 0 : static struct isl_tab *add_gbr_eq(struct isl_tab *tab, isl_int *eq)
3062 : {
3063 0 : if (!tab)
3064 0 : return NULL;
3065 :
3066 0 : if (isl_tab_extend_cons(tab, 2) < 0)
3067 0 : goto error;
3068 :
3069 0 : if (isl_tab_add_eq(tab, eq) < 0)
3070 0 : goto error;
3071 :
3072 0 : return tab;
3073 : error:
3074 0 : isl_tab_free(tab);
3075 0 : return NULL;
3076 : }
3077 :
3078 : /* Add the equality described by "eq" to the context.
3079 : * If "check" is set, then we check if the context is empty after
3080 : * adding the equality.
3081 : * If "update" is set, then we check if the samples are still valid.
3082 : *
3083 : * We do not explicitly add shifted copies of the equality to
3084 : * cgbr->shifted since they would conflict with each other.
3085 : * Instead, we directly mark cgbr->shifted empty.
3086 : */
3087 0 : static void context_gbr_add_eq(struct isl_context *context, isl_int *eq,
3088 : int check, int update)
3089 : {
3090 0 : struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3091 :
3092 0 : cgbr->tab = add_gbr_eq(cgbr->tab, eq);
3093 :
3094 0 : if (cgbr->shifted && !cgbr->shifted->empty && use_shifted(cgbr)) {
3095 0 : if (isl_tab_mark_empty(cgbr->shifted) < 0)
3096 0 : goto error;
3097 : }
3098 :
3099 0 : if (cgbr->cone && cgbr->cone->n_col != cgbr->cone->n_dead) {
3100 0 : if (isl_tab_extend_cons(cgbr->cone, 2) < 0)
3101 0 : goto error;
3102 0 : if (isl_tab_add_eq(cgbr->cone, eq) < 0)
3103 0 : goto error;
3104 : }
3105 :
3106 0 : if (check) {
3107 0 : int v = tab_has_valid_sample(cgbr->tab, eq, 1);
3108 0 : if (v < 0)
3109 0 : goto error;
3110 0 : if (!v)
3111 0 : check_gbr_integer_feasible(cgbr);
3112 : }
3113 0 : if (update)
3114 0 : cgbr->tab = check_samples(cgbr->tab, eq, 1);
3115 0 : return;
3116 : error:
3117 0 : isl_tab_free(cgbr->tab);
3118 0 : cgbr->tab = NULL;
3119 : }
3120 :
3121 0 : static void add_gbr_ineq(struct isl_context_gbr *cgbr, isl_int *ineq)
3122 : {
3123 0 : if (!cgbr->tab)
3124 0 : return;
3125 :
3126 0 : if (isl_tab_extend_cons(cgbr->tab, 1) < 0)
3127 0 : goto error;
3128 :
3129 0 : if (isl_tab_add_ineq(cgbr->tab, ineq) < 0)
3130 0 : goto error;
3131 :
3132 0 : if (cgbr->shifted && !cgbr->shifted->empty && use_shifted(cgbr)) {
3133 : int i;
3134 : unsigned dim;
3135 0 : dim = isl_basic_map_total_dim(cgbr->tab->bmap);
3136 :
3137 0 : if (isl_tab_extend_cons(cgbr->shifted, 1) < 0)
3138 0 : goto error;
3139 :
3140 0 : for (i = 0; i < dim; ++i) {
3141 0 : if (!isl_int_is_neg(ineq[1 + i]))
3142 0 : continue;
3143 0 : isl_int_add(ineq[0], ineq[0], ineq[1 + i]);
3144 : }
3145 :
3146 0 : if (isl_tab_add_ineq(cgbr->shifted, ineq) < 0)
3147 0 : goto error;
3148 :
3149 0 : for (i = 0; i < dim; ++i) {
3150 0 : if (!isl_int_is_neg(ineq[1 + i]))
3151 0 : continue;
3152 0 : isl_int_sub(ineq[0], ineq[0], ineq[1 + i]);
3153 : }
3154 : }
3155 :
3156 0 : if (cgbr->cone && cgbr->cone->n_col != cgbr->cone->n_dead) {
3157 0 : if (isl_tab_extend_cons(cgbr->cone, 1) < 0)
3158 0 : goto error;
3159 0 : if (isl_tab_add_ineq(cgbr->cone, ineq) < 0)
3160 0 : goto error;
3161 : }
3162 :
3163 0 : return;
3164 : error:
3165 0 : isl_tab_free(cgbr->tab);
3166 0 : cgbr->tab = NULL;
3167 : }
3168 :
3169 0 : static void context_gbr_add_ineq(struct isl_context *context, isl_int *ineq,
3170 : int check, int update)
3171 : {
3172 0 : struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3173 :
3174 0 : add_gbr_ineq(cgbr, ineq);
3175 0 : if (!cgbr->tab)
3176 0 : return;
3177 :
3178 0 : if (check) {
3179 0 : int v = tab_has_valid_sample(cgbr->tab, ineq, 0);
3180 0 : if (v < 0)
3181 0 : goto error;
3182 0 : if (!v)
3183 0 : check_gbr_integer_feasible(cgbr);
3184 : }
3185 0 : if (update)
3186 0 : cgbr->tab = check_samples(cgbr->tab, ineq, 0);
3187 0 : return;
3188 : error:
3189 0 : isl_tab_free(cgbr->tab);
3190 0 : cgbr->tab = NULL;
3191 : }
3192 :
3193 0 : static isl_stat context_gbr_add_ineq_wrap(void *user, isl_int *ineq)
3194 : {
3195 0 : struct isl_context *context = (struct isl_context *)user;
3196 0 : context_gbr_add_ineq(context, ineq, 0, 0);
3197 0 : return context->op->is_ok(context) ? isl_stat_ok : isl_stat_error;
3198 : }
3199 :
3200 0 : static enum isl_tab_row_sign context_gbr_ineq_sign(struct isl_context *context,
3201 : isl_int *ineq, int strict)
3202 : {
3203 0 : struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3204 0 : return tab_ineq_sign(cgbr->tab, ineq, strict);
3205 : }
3206 :
3207 : /* Check whether "ineq" can be added to the tableau without rendering
3208 : * it infeasible.
3209 : */
3210 0 : static int context_gbr_test_ineq(struct isl_context *context, isl_int *ineq)
3211 : {
3212 0 : struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3213 : struct isl_tab_undo *snap;
3214 0 : struct isl_tab_undo *shifted_snap = NULL;
3215 0 : struct isl_tab_undo *cone_snap = NULL;
3216 : int feasible;
3217 :
3218 0 : if (!cgbr->tab)
3219 0 : return -1;
3220 :
3221 0 : if (isl_tab_extend_cons(cgbr->tab, 1) < 0)
3222 0 : return -1;
3223 :
3224 0 : snap = isl_tab_snap(cgbr->tab);
3225 0 : if (cgbr->shifted)
3226 0 : shifted_snap = isl_tab_snap(cgbr->shifted);
3227 0 : if (cgbr->cone)
3228 0 : cone_snap = isl_tab_snap(cgbr->cone);
3229 0 : add_gbr_ineq(cgbr, ineq);
3230 0 : check_gbr_integer_feasible(cgbr);
3231 0 : if (!cgbr->tab)
3232 0 : return -1;
3233 0 : feasible = !cgbr->tab->empty;
3234 0 : if (isl_tab_rollback(cgbr->tab, snap) < 0)
3235 0 : return -1;
3236 0 : if (shifted_snap) {
3237 0 : if (isl_tab_rollback(cgbr->shifted, shifted_snap))
3238 0 : return -1;
3239 0 : } else if (cgbr->shifted) {
3240 0 : isl_tab_free(cgbr->shifted);
3241 0 : cgbr->shifted = NULL;
3242 : }
3243 0 : if (cone_snap) {
3244 0 : if (isl_tab_rollback(cgbr->cone, cone_snap))
3245 0 : return -1;
3246 0 : } else if (cgbr->cone) {
3247 0 : isl_tab_free(cgbr->cone);
3248 0 : cgbr->cone = NULL;
3249 : }
3250 :
3251 0 : return feasible;
3252 : }
3253 :
3254 : /* Return the column of the last of the variables associated to
3255 : * a column that has a non-zero coefficient.
3256 : * This function is called in a context where only coefficients
3257 : * of parameters or divs can be non-zero.
3258 : */
3259 0 : static int last_non_zero_var_col(struct isl_tab *tab, isl_int *p)
3260 : {
3261 : int i;
3262 : int col;
3263 :
3264 0 : if (tab->n_var == 0)
3265 0 : return -1;
3266 :
3267 0 : for (i = tab->n_var - 1; i >= 0; --i) {
3268 0 : if (i >= tab->n_param && i < tab->n_var - tab->n_div)
3269 0 : continue;
3270 0 : if (tab->var[i].is_row)
3271 0 : continue;
3272 0 : col = tab->var[i].index;
3273 0 : if (!isl_int_is_zero(p[col]))
3274 0 : return col;
3275 : }
3276 :
3277 0 : return -1;
3278 : }
3279 :
3280 : /* Look through all the recently added equalities in the context
3281 : * to see if we can propagate any of them to the main tableau.
3282 : *
3283 : * The newly added equalities in the context are encoded as pairs
3284 : * of inequalities starting at inequality "first".
3285 : *
3286 : * We tentatively add each of these equalities to the main tableau
3287 : * and if this happens to result in a row with a final coefficient
3288 : * that is one or negative one, we use it to kill a column
3289 : * in the main tableau. Otherwise, we discard the tentatively
3290 : * added row.
3291 : * This tentative addition of equality constraints turns
3292 : * on the undo facility of the tableau. Turn it off again
3293 : * at the end, assuming it was turned off to begin with.
3294 : *
3295 : * Return 0 on success and -1 on failure.
3296 : */
3297 0 : static int propagate_equalities(struct isl_context_gbr *cgbr,
3298 : struct isl_tab *tab, unsigned first)
3299 : {
3300 : int i;
3301 0 : struct isl_vec *eq = NULL;
3302 : isl_bool needs_undo;
3303 :
3304 0 : needs_undo = isl_tab_need_undo(tab);
3305 0 : if (needs_undo < 0)
3306 0 : goto error;
3307 0 : eq = isl_vec_alloc(tab->mat->ctx, 1 + tab->n_var);
3308 0 : if (!eq)
3309 0 : goto error;
3310 :
3311 0 : if (isl_tab_extend_cons(tab, (cgbr->tab->bmap->n_ineq - first)/2) < 0)
3312 0 : goto error;
3313 :
3314 0 : isl_seq_clr(eq->el + 1 + tab->n_param,
3315 0 : tab->n_var - tab->n_param - tab->n_div);
3316 0 : for (i = first; i < cgbr->tab->bmap->n_ineq; i += 2) {
3317 : int j;
3318 : int r;
3319 : struct isl_tab_undo *snap;
3320 0 : snap = isl_tab_snap(tab);
3321 :
3322 0 : isl_seq_cpy(eq->el, cgbr->tab->bmap->ineq[i], 1 + tab->n_param);
3323 0 : isl_seq_cpy(eq->el + 1 + tab->n_var - tab->n_div,
3324 0 : cgbr->tab->bmap->ineq[i] + 1 + tab->n_param,
3325 : tab->n_div);
3326 :
3327 0 : r = isl_tab_add_row(tab, eq->el);
3328 0 : if (r < 0)
3329 0 : goto error;
3330 0 : r = tab->con[r].index;
3331 0 : j = last_non_zero_var_col(tab, tab->mat->row[r] + 2 + tab->M);
3332 0 : if (j < 0 || j < tab->n_dead ||
3333 0 : !isl_int_is_one(tab->mat->row[r][0]) ||
3334 0 : (!isl_int_is_one(tab->mat->row[r][2 + tab->M + j]) &&
3335 0 : !isl_int_is_negone(tab->mat->row[r][2 + tab->M + j]))) {
3336 0 : if (isl_tab_rollback(tab, snap) < 0)
3337 0 : goto error;
3338 0 : continue;
3339 : }
3340 0 : if (isl_tab_pivot(tab, r, j) < 0)
3341 0 : goto error;
3342 0 : if (isl_tab_kill_col(tab, j) < 0)
3343 0 : goto error;
3344 :
3345 0 : if (restore_lexmin(tab) < 0)
3346 0 : goto error;
3347 : }
3348 :
3349 0 : if (!needs_undo)
3350 0 : isl_tab_clear_undo(tab);
3351 0 : isl_vec_free(eq);
3352 :
3353 0 : return 0;
3354 : error:
3355 0 : isl_vec_free(eq);
3356 0 : isl_tab_free(cgbr->tab);
3357 0 : cgbr->tab = NULL;
3358 0 : return -1;
3359 : }
3360 :
3361 0 : static int context_gbr_detect_equalities(struct isl_context *context,
3362 : struct isl_tab *tab)
3363 : {
3364 0 : struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3365 : unsigned n_ineq;
3366 :
3367 0 : if (!cgbr->cone) {
3368 0 : struct isl_basic_set *bset = isl_tab_peek_bset(cgbr->tab);
3369 0 : cgbr->cone = isl_tab_from_recession_cone(bset, 0);
3370 0 : if (!cgbr->cone)
3371 0 : goto error;
3372 0 : if (isl_tab_track_bset(cgbr->cone,
3373 : isl_basic_set_copy(bset)) < 0)
3374 0 : goto error;
3375 : }
3376 0 : if (isl_tab_detect_implicit_equalities(cgbr->cone) < 0)
3377 0 : goto error;
3378 :
3379 0 : n_ineq = cgbr->tab->bmap->n_ineq;
3380 0 : cgbr->tab = isl_tab_detect_equalities(cgbr->tab, cgbr->cone);
3381 0 : if (!cgbr->tab)
3382 0 : return -1;
3383 0 : if (cgbr->tab->bmap->n_ineq > n_ineq &&
3384 0 : propagate_equalities(cgbr, tab, n_ineq) < 0)
3385 0 : return -1;
3386 :
3387 0 : return 0;
3388 : error:
3389 0 : isl_tab_free(cgbr->tab);
3390 0 : cgbr->tab = NULL;
3391 0 : return -1;
3392 : }
3393 :
3394 0 : static int context_gbr_get_div(struct isl_context *context, struct isl_tab *tab,
3395 : struct isl_vec *div)
3396 : {
3397 0 : return get_div(tab, context, div);
3398 : }
3399 :
3400 0 : static isl_bool context_gbr_insert_div(struct isl_context *context, int pos,
3401 : __isl_keep isl_vec *div)
3402 : {
3403 0 : struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3404 0 : if (cgbr->cone) {
3405 : int r, n_div, o_div;
3406 :
3407 0 : n_div = isl_basic_map_dim(cgbr->cone->bmap, isl_dim_div);
3408 0 : o_div = cgbr->cone->n_var - n_div;
3409 :
3410 0 : if (isl_tab_extend_cons(cgbr->cone, 3) < 0)
3411 0 : return isl_bool_error;
3412 0 : if (isl_tab_extend_vars(cgbr->cone, 1) < 0)
3413 0 : return isl_bool_error;
3414 0 : if ((r = isl_tab_insert_var(cgbr->cone, pos)) <0)
3415 0 : return isl_bool_error;
3416 :
3417 0 : cgbr->cone->bmap = isl_basic_map_insert_div(cgbr->cone->bmap,
3418 : r - o_div, div);
3419 0 : if (!cgbr->cone->bmap)
3420 0 : return isl_bool_error;
3421 0 : if (isl_tab_push_var(cgbr->cone, isl_tab_undo_bmap_div,
3422 0 : &cgbr->cone->var[r]) < 0)
3423 0 : return isl_bool_error;
3424 : }
3425 0 : return context_tab_insert_div(cgbr->tab, pos, div,
3426 : context_gbr_add_ineq_wrap, context);
3427 : }
3428 :
3429 0 : static int context_gbr_best_split(struct isl_context *context,
3430 : struct isl_tab *tab)
3431 : {
3432 0 : struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3433 : struct isl_tab_undo *snap;
3434 : int r;
3435 :
3436 0 : snap = isl_tab_snap(cgbr->tab);
3437 0 : r = best_split(tab, cgbr->tab);
3438 :
3439 0 : if (r >= 0 && isl_tab_rollback(cgbr->tab, snap) < 0)
3440 0 : return -1;
3441 :
3442 0 : return r;
3443 : }
3444 :
3445 0 : static int context_gbr_is_empty(struct isl_context *context)
3446 : {
3447 0 : struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3448 0 : if (!cgbr->tab)
3449 0 : return -1;
3450 0 : return cgbr->tab->empty;
3451 : }
3452 :
3453 : struct isl_gbr_tab_undo {
3454 : struct isl_tab_undo *tab_snap;
3455 : struct isl_tab_undo *shifted_snap;
3456 : struct isl_tab_undo *cone_snap;
3457 : };
3458 :
3459 0 : static void *context_gbr_save(struct isl_context *context)
3460 : {
3461 0 : struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3462 : struct isl_gbr_tab_undo *snap;
3463 :
3464 0 : if (!cgbr->tab)
3465 0 : return NULL;
3466 :
3467 0 : snap = isl_alloc_type(cgbr->tab->mat->ctx, struct isl_gbr_tab_undo);
3468 0 : if (!snap)
3469 0 : return NULL;
3470 :
3471 0 : snap->tab_snap = isl_tab_snap(cgbr->tab);
3472 0 : if (isl_tab_save_samples(cgbr->tab) < 0)
3473 0 : goto error;
3474 :
3475 0 : if (cgbr->shifted)
3476 0 : snap->shifted_snap = isl_tab_snap(cgbr->shifted);
3477 : else
3478 0 : snap->shifted_snap = NULL;
3479 :
3480 0 : if (cgbr->cone)
3481 0 : snap->cone_snap = isl_tab_snap(cgbr->cone);
3482 : else
3483 0 : snap->cone_snap = NULL;
3484 :
3485 0 : return snap;
3486 : error:
3487 0 : free(snap);
3488 0 : return NULL;
3489 : }
3490 :
3491 0 : static void context_gbr_restore(struct isl_context *context, void *save)
3492 : {
3493 0 : struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3494 0 : struct isl_gbr_tab_undo *snap = (struct isl_gbr_tab_undo *)save;
3495 0 : if (!snap)
3496 0 : goto error;
3497 0 : if (isl_tab_rollback(cgbr->tab, snap->tab_snap) < 0)
3498 0 : goto error;
3499 :
3500 0 : if (snap->shifted_snap) {
3501 0 : if (isl_tab_rollback(cgbr->shifted, snap->shifted_snap) < 0)
3502 0 : goto error;
3503 0 : } else if (cgbr->shifted) {
3504 0 : isl_tab_free(cgbr->shifted);
3505 0 : cgbr->shifted = NULL;
3506 : }
3507 :
3508 0 : if (snap->cone_snap) {
3509 0 : if (isl_tab_rollback(cgbr->cone, snap->cone_snap) < 0)
3510 0 : goto error;
3511 0 : } else if (cgbr->cone) {
3512 0 : isl_tab_free(cgbr->cone);
3513 0 : cgbr->cone = NULL;
3514 : }
3515 :
3516 0 : free(snap);
3517 :
3518 0 : return;
3519 : error:
3520 0 : free(snap);
3521 0 : isl_tab_free(cgbr->tab);
3522 0 : cgbr->tab = NULL;
3523 : }
3524 :
3525 0 : static void context_gbr_discard(void *save)
3526 : {
3527 0 : struct isl_gbr_tab_undo *snap = (struct isl_gbr_tab_undo *)save;
3528 0 : free(snap);
3529 0 : }
3530 :
3531 0 : static int context_gbr_is_ok(struct isl_context *context)
3532 : {
3533 0 : struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3534 0 : return !!cgbr->tab;
3535 : }
3536 :
3537 0 : static void context_gbr_invalidate(struct isl_context *context)
3538 : {
3539 0 : struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3540 0 : isl_tab_free(cgbr->tab);
3541 0 : cgbr->tab = NULL;
3542 0 : }
3543 :
3544 0 : static __isl_null struct isl_context *context_gbr_free(
3545 : struct isl_context *context)
3546 : {
3547 0 : struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3548 0 : isl_tab_free(cgbr->tab);
3549 0 : isl_tab_free(cgbr->shifted);
3550 0 : isl_tab_free(cgbr->cone);
3551 0 : free(cgbr);
3552 :
3553 0 : return NULL;
3554 : }
3555 :
3556 : struct isl_context_op isl_context_gbr_op = {
3557 : context_gbr_detect_nonnegative_parameters,
3558 : context_gbr_peek_basic_set,
3559 : context_gbr_peek_tab,
3560 : context_gbr_add_eq,
3561 : context_gbr_add_ineq,
3562 : context_gbr_ineq_sign,
3563 : context_gbr_test_ineq,
3564 : context_gbr_get_div,
3565 : context_gbr_insert_div,
3566 : context_gbr_detect_equalities,
3567 : context_gbr_best_split,
3568 : context_gbr_is_empty,
3569 : context_gbr_is_ok,
3570 : context_gbr_save,
3571 : context_gbr_restore,
3572 : context_gbr_discard,
3573 : context_gbr_invalidate,
3574 : context_gbr_free,
3575 : };
3576 :
3577 0 : static struct isl_context *isl_context_gbr_alloc(__isl_keep isl_basic_set *dom)
3578 : {
3579 : struct isl_context_gbr *cgbr;
3580 :
3581 0 : if (!dom)
3582 0 : return NULL;
3583 :
3584 0 : cgbr = isl_calloc_type(dom->ctx, struct isl_context_gbr);
3585 0 : if (!cgbr)
3586 0 : return NULL;
3587 :
3588 0 : cgbr->context.op = &isl_context_gbr_op;
3589 :
3590 0 : cgbr->shifted = NULL;
3591 0 : cgbr->cone = NULL;
3592 0 : cgbr->tab = isl_tab_from_basic_set(dom, 1);
3593 0 : cgbr->tab = isl_tab_init_samples(cgbr->tab);
3594 0 : if (!cgbr->tab)
3595 0 : goto error;
3596 0 : check_gbr_integer_feasible(cgbr);
3597 :
3598 0 : return &cgbr->context;
3599 : error:
3600 0 : cgbr->context.op->free(&cgbr->context);
3601 0 : return NULL;
3602 : }
3603 :
3604 : /* Allocate a context corresponding to "dom".
3605 : * The representation specific fields are initialized by
3606 : * isl_context_lex_alloc or isl_context_gbr_alloc.
3607 : * The shared "n_unknown" field is initialized to the number
3608 : * of final unknown integer divisions in "dom".
3609 : */
3610 0 : static struct isl_context *isl_context_alloc(__isl_keep isl_basic_set *dom)
3611 : {
3612 : struct isl_context *context;
3613 : int first;
3614 :
3615 0 : if (!dom)
3616 0 : return NULL;
3617 :
3618 0 : if (dom->ctx->opt->context == ISL_CONTEXT_LEXMIN)
3619 0 : context = isl_context_lex_alloc(dom);
3620 : else
3621 0 : context = isl_context_gbr_alloc(dom);
3622 :
3623 0 : if (!context)
3624 0 : return NULL;
3625 :
3626 0 : first = isl_basic_set_first_unknown_div(dom);
3627 0 : if (first < 0)
3628 0 : return context->op->free(context);
3629 0 : context->n_unknown = isl_basic_set_dim(dom, isl_dim_div) - first;
3630 :
3631 0 : return context;
3632 : }
3633 :
3634 : /* Initialize some common fields of "sol", which keeps track
3635 : * of the solution of an optimization problem on "bmap" over
3636 : * the domain "dom".
3637 : * If "max" is set, then a maximization problem is being solved, rather than
3638 : * a minimization problem, which means that the variables in the
3639 : * tableau have value "M - x" rather than "M + x".
3640 : */
3641 0 : static isl_stat sol_init(struct isl_sol *sol, __isl_keep isl_basic_map *bmap,
3642 : __isl_keep isl_basic_set *dom, int max)
3643 : {
3644 0 : sol->rational = ISL_F_ISSET(bmap, ISL_BASIC_MAP_RATIONAL);
3645 0 : sol->dec_level.callback.run = &sol_dec_level_wrap;
3646 0 : sol->dec_level.sol = sol;
3647 0 : sol->max = max;
3648 0 : sol->n_out = isl_basic_map_dim(bmap, isl_dim_out);
3649 0 : sol->space = isl_basic_map_get_space(bmap);
3650 :
3651 0 : sol->context = isl_context_alloc(dom);
3652 0 : if (!sol->space || !sol->context)
3653 0 : return isl_stat_error;
3654 :
3655 0 : return isl_stat_ok;
3656 : }
3657 :
3658 : /* Construct an isl_sol_map structure for accumulating the solution.
3659 : * If track_empty is set, then we also keep track of the parts
3660 : * of the context where there is no solution.
3661 : * If max is set, then we are solving a maximization, rather than
3662 : * a minimization problem, which means that the variables in the
3663 : * tableau have value "M - x" rather than "M + x".
3664 : */
3665 0 : static struct isl_sol *sol_map_init(__isl_keep isl_basic_map *bmap,
3666 : __isl_take isl_basic_set *dom, int track_empty, int max)
3667 : {
3668 0 : struct isl_sol_map *sol_map = NULL;
3669 : isl_space *space;
3670 :
3671 0 : if (!bmap)
3672 0 : goto error;
3673 :
3674 0 : sol_map = isl_calloc_type(bmap->ctx, struct isl_sol_map);
3675 0 : if (!sol_map)
3676 0 : goto error;
3677 :
3678 0 : sol_map->sol.free = &sol_map_free;
3679 0 : if (sol_init(&sol_map->sol, bmap, dom, max) < 0)
3680 0 : goto error;
3681 0 : sol_map->sol.add = &sol_map_add_wrap;
3682 0 : sol_map->sol.add_empty = track_empty ? &sol_map_add_empty_wrap : NULL;
3683 0 : space = isl_space_copy(sol_map->sol.space);
3684 0 : sol_map->map = isl_map_alloc_space(space, 1, ISL_MAP_DISJOINT);
3685 0 : if (!sol_map->map)
3686 0 : goto error;
3687 :
3688 0 : if (track_empty) {
3689 0 : sol_map->empty = isl_set_alloc_space(isl_basic_set_get_space(dom),
3690 : 1, ISL_SET_DISJOINT);
3691 0 : if (!sol_map->empty)
3692 0 : goto error;
3693 : }
3694 :
3695 0 : isl_basic_set_free(dom);
3696 0 : return &sol_map->sol;
3697 : error:
3698 0 : isl_basic_set_free(dom);
3699 0 : sol_free(&sol_map->sol);
3700 0 : return NULL;
3701 : }
3702 :
3703 : /* Check whether all coefficients of (non-parameter) variables
3704 : * are non-positive, meaning that no pivots can be performed on the row.
3705 : */
3706 0 : static int is_critical(struct isl_tab *tab, int row)
3707 : {
3708 : int j;
3709 0 : unsigned off = 2 + tab->M;
3710 :
3711 0 : for (j = tab->n_dead; j < tab->n_col; ++j) {
3712 0 : if (col_is_parameter_var(tab, j))
3713 0 : continue;
3714 :
3715 0 : if (isl_int_is_pos(tab->mat->row[row][off + j]))
3716 0 : return 0;
3717 : }
3718 :
3719 0 : return 1;
3720 : }
3721 :
3722 : /* Check whether the inequality represented by vec is strict over the integers,
3723 : * i.e., there are no integer values satisfying the constraint with
3724 : * equality. This happens if the gcd of the coefficients is not a divisor
3725 : * of the constant term. If so, scale the constraint down by the gcd
3726 : * of the coefficients.
3727 : */
3728 0 : static int is_strict(struct isl_vec *vec)
3729 : {
3730 : isl_int gcd;
3731 0 : int strict = 0;
3732 :
3733 0 : isl_int_init(gcd);
3734 0 : isl_seq_gcd(vec->el + 1, vec->size - 1, &gcd);
3735 0 : if (!isl_int_is_one(gcd)) {
3736 0 : strict = !isl_int_is_divisible_by(vec->el[0], gcd);
3737 0 : isl_int_fdiv_q(vec->el[0], vec->el[0], gcd);
3738 0 : isl_seq_scale_down(vec->el + 1, vec->el + 1, gcd, vec->size-1);
3739 : }
3740 0 : isl_int_clear(gcd);
3741 :
3742 0 : return strict;
3743 : }
3744 :
3745 : /* Determine the sign of the given row of the main tableau.
3746 : * The result is one of
3747 : * isl_tab_row_pos: always non-negative; no pivot needed
3748 : * isl_tab_row_neg: always non-positive; pivot
3749 : * isl_tab_row_any: can be both positive and negative; split
3750 : *
3751 : * We first handle some simple cases
3752 : * - the row sign may be known already
3753 : * - the row may be obviously non-negative
3754 : * - the parametric constant may be equal to that of another row
3755 : * for which we know the sign. This sign will be either "pos" or
3756 : * "any". If it had been "neg" then we would have pivoted before.
3757 : *
3758 : * If none of these cases hold, we check the value of the row for each
3759 : * of the currently active samples. Based on the signs of these values
3760 : * we make an initial determination of the sign of the row.
3761 : *
3762 : * all zero -> unk(nown)
3763 : * all non-negative -> pos
3764 : * all non-positive -> neg
3765 : * both negative and positive -> all
3766 : *
3767 : * If we end up with "all", we are done.
3768 : * Otherwise, we perform a check for positive and/or negative
3769 : * values as follows.
3770 : *
3771 : * samples neg unk pos
3772 : * <0 ? Y N Y N
3773 : * pos any pos
3774 : * >0 ? Y N Y N
3775 : * any neg any neg
3776 : *
3777 : * There is no special sign for "zero", because we can usually treat zero
3778 : * as either non-negative or non-positive, whatever works out best.
3779 : * However, if the row is "critical", meaning that pivoting is impossible
3780 : * then we don't want to limp zero with the non-positive case, because
3781 : * then we we would lose the solution for those values of the parameters
3782 : * where the value of the row is zero. Instead, we treat 0 as non-negative
3783 : * ensuring a split if the row can attain both zero and negative values.
3784 : * The same happens when the original constraint was one that could not
3785 : * be satisfied with equality by any integer values of the parameters.
3786 : * In this case, we normalize the constraint, but then a value of zero
3787 : * for the normalized constraint is actually a positive value for the
3788 : * original constraint, so again we need to treat zero as non-negative.
3789 : * In both these cases, we have the following decision tree instead:
3790 : *
3791 : * all non-negative -> pos
3792 : * all negative -> neg
3793 : * both negative and non-negative -> all
3794 : *
3795 : * samples neg pos
3796 : * <0 ? Y N
3797 : * any pos
3798 : * >=0 ? Y N
3799 : * any neg
3800 : */
3801 0 : static enum isl_tab_row_sign row_sign(struct isl_tab *tab,
3802 : struct isl_sol *sol, int row)
3803 : {
3804 0 : struct isl_vec *ineq = NULL;
3805 0 : enum isl_tab_row_sign res = isl_tab_row_unknown;
3806 : int critical;
3807 : int strict;
3808 : int row2;
3809 :
3810 0 : if (tab->row_sign[row] != isl_tab_row_unknown)
3811 0 : return tab->row_sign[row];
3812 0 : if (is_obviously_nonneg(tab, row))
3813 0 : return isl_tab_row_pos;
3814 0 : for (row2 = tab->n_redundant; row2 < tab->n_row; ++row2) {
3815 0 : if (tab->row_sign[row2] == isl_tab_row_unknown)
3816 0 : continue;
3817 0 : if (identical_parameter_line(tab, row, row2))
3818 0 : return tab->row_sign[row2];
3819 : }
3820 :
3821 0 : critical = is_critical(tab, row);
3822 :
3823 0 : ineq = get_row_parameter_ineq(tab, row);
3824 0 : if (!ineq)
3825 0 : goto error;
3826 :
3827 0 : strict = is_strict(ineq);
3828 :
3829 0 : res = sol->context->op->ineq_sign(sol->context, ineq->el,
3830 : critical || strict);
3831 :
3832 0 : if (res == isl_tab_row_unknown || res == isl_tab_row_pos) {
3833 : /* test for negative values */
3834 : int feasible;
3835 0 : isl_seq_neg(ineq->el, ineq->el, ineq->size);
3836 0 : isl_int_sub_ui(ineq->el[0], ineq->el[0], 1);
3837 :
3838 0 : feasible = sol->context->op->test_ineq(sol->context, ineq->el);
3839 0 : if (feasible < 0)
3840 0 : goto error;
3841 0 : if (!feasible)
3842 0 : res = isl_tab_row_pos;
3843 : else
3844 0 : res = (res == isl_tab_row_unknown) ? isl_tab_row_neg
3845 : : isl_tab_row_any;
3846 0 : if (res == isl_tab_row_neg) {
3847 0 : isl_seq_neg(ineq->el, ineq->el, ineq->size);
3848 0 : isl_int_sub_ui(ineq->el[0], ineq->el[0], 1);
3849 : }
3850 : }
3851 :
3852 0 : if (res == isl_tab_row_neg) {
3853 : /* test for positive values */
3854 : int feasible;
3855 0 : if (!critical && !strict)
3856 0 : isl_int_sub_ui(ineq->el[0], ineq->el[0], 1);
3857 :
3858 0 : feasible = sol->context->op->test_ineq(sol->context, ineq->el);
3859 0 : if (feasible < 0)
3860 0 : goto error;
3861 0 : if (feasible)
3862 0 : res = isl_tab_row_any;
3863 : }
3864 :
3865 0 : isl_vec_free(ineq);
3866 0 : return res;
3867 : error:
3868 0 : isl_vec_free(ineq);
3869 0 : return isl_tab_row_unknown;
3870 : }
3871 :
3872 : static void find_solutions(struct isl_sol *sol, struct isl_tab *tab);
3873 :
3874 : /* Find solutions for values of the parameters that satisfy the given
3875 : * inequality.
3876 : *
3877 : * We currently take a snapshot of the context tableau that is reset
3878 : * when we return from this function, while we make a copy of the main
3879 : * tableau, leaving the original main tableau untouched.
3880 : * These are fairly arbitrary choices. Making a copy also of the context
3881 : * tableau would obviate the need to undo any changes made to it later,
3882 : * while taking a snapshot of the main tableau could reduce memory usage.
3883 : * If we were to switch to taking a snapshot of the main tableau,
3884 : * we would have to keep in mind that we need to save the row signs
3885 : * and that we need to do this before saving the current basis
3886 : * such that the basis has been restore before we restore the row signs.
3887 : */
3888 0 : static void find_in_pos(struct isl_sol *sol, struct isl_tab *tab, isl_int *ineq)
3889 : {
3890 : void *saved;
3891 :
3892 0 : if (!sol->context)
3893 0 : goto error;
3894 0 : saved = sol->context->op->save(sol->context);
3895 :
3896 0 : tab = isl_tab_dup(tab);
3897 0 : if (!tab)
3898 0 : goto error;
3899 :
3900 0 : sol->context->op->add_ineq(sol->context, ineq, 0, 1);
3901 :
3902 0 : find_solutions(sol, tab);
3903 :
3904 0 : if (!sol->error)
3905 0 : sol->context->op->restore(sol->context, saved);
3906 : else
3907 0 : sol->context->op->discard(saved);
3908 0 : return;
3909 : error:
3910 0 : sol->error = 1;
3911 : }
3912 :
3913 : /* Record the absence of solutions for those values of the parameters
3914 : * that do not satisfy the given inequality with equality.
3915 : */
3916 0 : static void no_sol_in_strict(struct isl_sol *sol,
3917 : struct isl_tab *tab, struct isl_vec *ineq)
3918 : {
3919 : int empty;
3920 : void *saved;
3921 :
3922 0 : if (!sol->context || sol->error)
3923 : goto error;
3924 0 : saved = sol->context->op->save(sol->context);
3925 :
3926 0 : isl_int_sub_ui(ineq->el[0], ineq->el[0], 1);
3927 :
3928 0 : sol->context->op->add_ineq(sol->context, ineq->el, 1, 0);
3929 0 : if (!sol->context)
3930 0 : goto error;
3931 :
3932 0 : empty = tab->empty;
3933 0 : tab->empty = 1;
3934 0 : sol_add(sol, tab);
3935 0 : tab->empty = empty;
3936 :
3937 0 : isl_int_add_ui(ineq->el[0], ineq->el[0], 1);
3938 :
3939 0 : sol->context->op->restore(sol->context, saved);
3940 0 : return;
3941 : error:
3942 0 : sol->error = 1;
3943 : }
3944 :
3945 : /* Reset all row variables that are marked to have a sign that may
3946 : * be both positive and negative to have an unknown sign.
3947 : */
3948 0 : static void reset_any_to_unknown(struct isl_tab *tab)
3949 : {
3950 : int row;
3951 :
3952 0 : for (row = tab->n_redundant; row < tab->n_row; ++row) {
3953 0 : if (!isl_tab_var_from_row(tab, row)->is_nonneg)
3954 0 : continue;
3955 0 : if (tab->row_sign[row] == isl_tab_row_any)
3956 0 : tab->row_sign[row] = isl_tab_row_unknown;
3957 : }
3958 0 : }
3959 :
3960 : /* Compute the lexicographic minimum of the set represented by the main
3961 : * tableau "tab" within the context "sol->context_tab".
3962 : * On entry the sample value of the main tableau is lexicographically
3963 : * less than or equal to this lexicographic minimum.
3964 : * Pivots are performed until a feasible point is found, which is then
3965 : * necessarily equal to the minimum, or until the tableau is found to
3966 : * be infeasible. Some pivots may need to be performed for only some
3967 : * feasible values of the context tableau. If so, the context tableau
3968 : * is split into a part where the pivot is needed and a part where it is not.
3969 : *
3970 : * Whenever we enter the main loop, the main tableau is such that no
3971 : * "obvious" pivots need to be performed on it, where "obvious" means
3972 : * that the given row can be seen to be negative without looking at
3973 : * the context tableau. In particular, for non-parametric problems,
3974 : * no pivots need to be performed on the main tableau.
3975 : * The caller of find_solutions is responsible for making this property
3976 : * hold prior to the first iteration of the loop, while restore_lexmin
3977 : * is called before every other iteration.
3978 : *
3979 : * Inside the main loop, we first examine the signs of the rows of
3980 : * the main tableau within the context of the context tableau.
3981 : * If we find a row that is always non-positive for all values of
3982 : * the parameters satisfying the context tableau and negative for at
3983 : * least one value of the parameters, we perform the appropriate pivot
3984 : * and start over. An exception is the case where no pivot can be
3985 : * performed on the row. In this case, we require that the sign of
3986 : * the row is negative for all values of the parameters (rather than just
3987 : * non-positive). This special case is handled inside row_sign, which
3988 : * will say that the row can have any sign if it determines that it can
3989 : * attain both negative and zero values.
3990 : *
3991 : * If we can't find a row that always requires a pivot, but we can find
3992 : * one or more rows that require a pivot for some values of the parameters
3993 : * (i.e., the row can attain both positive and negative signs), then we split
3994 : * the context tableau into two parts, one where we force the sign to be
3995 : * non-negative and one where we force is to be negative.
3996 : * The non-negative part is handled by a recursive call (through find_in_pos).
3997 : * Upon returning from this call, we continue with the negative part and
3998 : * perform the required pivot.
3999 : *
4000 : * If no such rows can be found, all rows are non-negative and we have
4001 : * found a (rational) feasible point. If we only wanted a rational point
4002 : * then we are done.
4003 : * Otherwise, we check if all values of the sample point of the tableau
4004 : * are integral for the variables. If so, we have found the minimal
4005 : * integral point and we are done.
4006 : * If the sample point is not integral, then we need to make a distinction
4007 : * based on whether the constant term is non-integral or the coefficients
4008 : * of the parameters. Furthermore, in order to decide how to handle
4009 : * the non-integrality, we also need to know whether the coefficients
4010 : * of the other columns in the tableau are integral. This leads
4011 : * to the following table. The first two rows do not correspond
4012 : * to a non-integral sample point and are only mentioned for completeness.
4013 : *
4014 : * constant parameters other
4015 : *
4016 : * int int int |
4017 : * int int rat | -> no problem
4018 : *
4019 : * rat int int -> fail
4020 : *
4021 : * rat int rat -> cut
4022 : *
4023 : * int rat rat |
4024 : * rat rat rat | -> parametric cut
4025 : *
4026 : * int rat int |
4027 : * rat rat int | -> split context
4028 : *
4029 : * If the parametric constant is completely integral, then there is nothing
4030 : * to be done. If the constant term is non-integral, but all the other
4031 : * coefficient are integral, then there is nothing that can be done
4032 : * and the tableau has no integral solution.
4033 : * If, on the other hand, one or more of the other columns have rational
4034 : * coefficients, but the parameter coefficients are all integral, then
4035 : * we can perform a regular (non-parametric) cut.
4036 : * Finally, if there is any parameter coefficient that is non-integral,
4037 : * then we need to involve the context tableau. There are two cases here.
4038 : * If at least one other column has a rational coefficient, then we
4039 : * can perform a parametric cut in the main tableau by adding a new
4040 : * integer division in the context tableau.
4041 : * If all other columns have integral coefficients, then we need to
4042 : * enforce that the rational combination of parameters (c + \sum a_i y_i)/m
4043 : * is always integral. We do this by introducing an integer division
4044 : * q = floor((c + \sum a_i y_i)/m) and stipulating that its argument should
4045 : * always be integral in the context tableau, i.e., m q = c + \sum a_i y_i.
4046 : * Since q is expressed in the tableau as
4047 : * c + \sum a_i y_i - m q >= 0
4048 : * -c - \sum a_i y_i + m q + m - 1 >= 0
4049 : * it is sufficient to add the inequality
4050 : * -c - \sum a_i y_i + m q >= 0
4051 : * In the part of the context where this inequality does not hold, the
4052 : * main tableau is marked as being empty.
4053 : */
4054 0 : static void find_solutions(struct isl_sol *sol, struct isl_tab *tab)
4055 : {
4056 : struct isl_context *context;
4057 : int r;
4058 :
4059 0 : if (!tab || sol->error)
4060 : goto error;
4061 :
4062 0 : context = sol->context;
4063 :
4064 0 : if (tab->empty)
4065 0 : goto done;
4066 0 : if (context->op->is_empty(context))
4067 0 : goto done;
4068 :
4069 0 : for (r = 0; r >= 0 && tab && !tab->empty; r = restore_lexmin(tab)) {
4070 : int flags;
4071 : int row;
4072 : enum isl_tab_row_sign sgn;
4073 0 : int split = -1;
4074 0 : int n_split = 0;
4075 :
4076 0 : for (row = tab->n_redundant; row < tab->n_row; ++row) {
4077 0 : if (!isl_tab_var_from_row(tab, row)->is_nonneg)
4078 0 : continue;
4079 0 : sgn = row_sign(tab, sol, row);
4080 0 : if (!sgn)
4081 0 : goto error;
4082 0 : tab->row_sign[row] = sgn;
4083 0 : if (sgn == isl_tab_row_any)
4084 0 : n_split++;
4085 0 : if (sgn == isl_tab_row_any && split == -1)
4086 0 : split = row;
4087 0 : if (sgn == isl_tab_row_neg)
4088 0 : break;
4089 : }
4090 0 : if (row < tab->n_row)
4091 0 : continue;
4092 0 : if (split != -1) {
4093 : struct isl_vec *ineq;
4094 0 : if (n_split != 1)
4095 0 : split = context->op->best_split(context, tab);
4096 0 : if (split < 0)
4097 0 : goto error;
4098 0 : ineq = get_row_parameter_ineq(tab, split);
4099 0 : if (!ineq)
4100 0 : goto error;
4101 0 : is_strict(ineq);
4102 0 : reset_any_to_unknown(tab);
4103 0 : tab->row_sign[split] = isl_tab_row_pos;
4104 0 : sol_inc_level(sol);
4105 0 : find_in_pos(sol, tab, ineq->el);
4106 0 : tab->row_sign[split] = isl_tab_row_neg;
4107 0 : isl_seq_neg(ineq->el, ineq->el, ineq->size);
4108 0 : isl_int_sub_ui(ineq->el[0], ineq->el[0], 1);
4109 0 : if (!sol->error)
4110 0 : context->op->add_ineq(context, ineq->el, 0, 1);
4111 0 : isl_vec_free(ineq);
4112 0 : if (sol->error)
4113 0 : goto error;
4114 0 : continue;
4115 : }
4116 0 : if (tab->rational)
4117 0 : break;
4118 0 : row = first_non_integer_row(tab, &flags);
4119 0 : if (row < 0)
4120 0 : break;
4121 0 : if (ISL_FL_ISSET(flags, I_PAR)) {
4122 0 : if (ISL_FL_ISSET(flags, I_VAR)) {
4123 0 : if (isl_tab_mark_empty(tab) < 0)
4124 0 : goto error;
4125 0 : break;
4126 : }
4127 0 : row = add_cut(tab, row);
4128 0 : } else if (ISL_FL_ISSET(flags, I_VAR)) {
4129 : struct isl_vec *div;
4130 : struct isl_vec *ineq;
4131 : int d;
4132 0 : div = get_row_split_div(tab, row);
4133 0 : if (!div)
4134 0 : goto error;
4135 0 : d = context->op->get_div(context, tab, div);
4136 0 : isl_vec_free(div);
4137 0 : if (d < 0)
4138 0 : goto error;
4139 0 : ineq = ineq_for_div(context->op->peek_basic_set(context), d);
4140 0 : if (!ineq)
4141 0 : goto error;
4142 0 : sol_inc_level(sol);
4143 0 : no_sol_in_strict(sol, tab, ineq);
4144 0 : isl_seq_neg(ineq->el, ineq->el, ineq->size);
4145 0 : context->op->add_ineq(context, ineq->el, 1, 1);
4146 0 : isl_vec_free(ineq);
4147 0 : if (sol->error || !context->op->is_ok(context))
4148 : goto error;
4149 0 : tab = set_row_cst_to_div(tab, row, d);
4150 0 : if (context->op->is_empty(context))
4151 0 : break;
4152 : } else
4153 0 : row = add_parametric_cut(tab, row, context);
4154 0 : if (row < 0)
4155 0 : goto error;
4156 : }
4157 0 : if (r < 0)
4158 0 : goto error;
4159 : done:
4160 0 : sol_add(sol, tab);
4161 0 : isl_tab_free(tab);
4162 0 : return;
4163 : error:
4164 0 : isl_tab_free(tab);
4165 0 : sol->error = 1;
4166 : }
4167 :
4168 : /* Does "sol" contain a pair of partial solutions that could potentially
4169 : * be merged?
4170 : *
4171 : * We currently only check that "sol" is not in an error state
4172 : * and that there are at least two partial solutions of which the final two
4173 : * are defined at the same level.
4174 : */
4175 0 : static int sol_has_mergeable_solutions(struct isl_sol *sol)
4176 : {
4177 0 : if (sol->error)
4178 0 : return 0;
4179 0 : if (!sol->partial)
4180 0 : return 0;
4181 0 : if (!sol->partial->next)
4182 0 : return 0;
4183 0 : return sol->partial->level == sol->partial->next->level;
4184 : }
4185 :
4186 : /* Compute the lexicographic minimum of the set represented by the main
4187 : * tableau "tab" within the context "sol->context_tab".
4188 : *
4189 : * As a preprocessing step, we first transfer all the purely parametric
4190 : * equalities from the main tableau to the context tableau, i.e.,
4191 : * parameters that have been pivoted to a row.
4192 : * These equalities are ignored by the main algorithm, because the
4193 : * corresponding rows may not be marked as being non-negative.
4194 : * In parts of the context where the added equality does not hold,
4195 : * the main tableau is marked as being empty.
4196 : *
4197 : * Before we embark on the actual computation, we save a copy
4198 : * of the context. When we return, we check if there are any
4199 : * partial solutions that can potentially be merged. If so,
4200 : * we perform a rollback to the initial state of the context.
4201 : * The merging of partial solutions happens inside calls to
4202 : * sol_dec_level that are pushed onto the undo stack of the context.
4203 : * If there are no partial solutions that can potentially be merged
4204 : * then the rollback is skipped as it would just be wasted effort.
4205 : */
4206 0 : static void find_solutions_main(struct isl_sol *sol, struct isl_tab *tab)
4207 : {
4208 : int row;
4209 : void *saved;
4210 :
4211 0 : if (!tab)
4212 0 : goto error;
4213 :
4214 0 : sol->level = 0;
4215 :
4216 0 : for (row = tab->n_redundant; row < tab->n_row; ++row) {
4217 : int p;
4218 : struct isl_vec *eq;
4219 :
4220 0 : if (!row_is_parameter_var(tab, row))
4221 0 : continue;
4222 0 : if (tab->row_var[row] < tab->n_param)
4223 0 : p = tab->row_var[row];
4224 : else
4225 0 : p = tab->row_var[row]
4226 0 : + tab->n_param - (tab->n_var - tab->n_div);
4227 :
4228 0 : eq = isl_vec_alloc(tab->mat->ctx, 1+tab->n_param+tab->n_div);
4229 0 : if (!eq)
4230 0 : goto error;
4231 0 : get_row_parameter_line(tab, row, eq->el);
4232 0 : isl_int_neg(eq->el[1 + p], tab->mat->row[row][0]);
4233 0 : eq = isl_vec_normalize(eq);
4234 :
4235 0 : sol_inc_level(sol);
4236 0 : no_sol_in_strict(sol, tab, eq);
4237 :
4238 0 : isl_seq_neg(eq->el, eq->el, eq->size);
4239 0 : sol_inc_level(sol);
4240 0 : no_sol_in_strict(sol, tab, eq);
4241 0 : isl_seq_neg(eq->el, eq->el, eq->size);
4242 :
4243 0 : sol->context->op->add_eq(sol->context, eq->el, 1, 1);
4244 :
4245 0 : isl_vec_free(eq);
4246 :
4247 0 : if (isl_tab_mark_redundant(tab, row) < 0)
4248 0 : goto error;
4249 :
4250 0 : if (sol->context->op->is_empty(sol->context))
4251 0 : break;
4252 :
4253 0 : row = tab->n_redundant - 1;
4254 : }
4255 :
4256 0 : saved = sol->context->op->save(sol->context);
4257 :
4258 0 : find_solutions(sol, tab);
4259 :
4260 0 : if (sol_has_mergeable_solutions(sol))
4261 0 : sol->context->op->restore(sol->context, saved);
4262 : else
4263 0 : sol->context->op->discard(saved);
4264 :
4265 0 : sol->level = 0;
4266 0 : sol_pop(sol);
4267 :
4268 0 : return;
4269 : error:
4270 0 : isl_tab_free(tab);
4271 0 : sol->error = 1;
4272 : }
4273 :
4274 : /* Check if integer division "div" of "dom" also occurs in "bmap".
4275 : * If so, return its position within the divs.
4276 : * If not, return -1.
4277 : */
4278 0 : static int find_context_div(struct isl_basic_map *bmap,
4279 : struct isl_basic_set *dom, unsigned div)
4280 : {
4281 : int i;
4282 0 : unsigned b_dim = isl_space_dim(bmap->dim, isl_dim_all);
4283 0 : unsigned d_dim = isl_space_dim(dom->dim, isl_dim_all);
4284 :
4285 0 : if (isl_int_is_zero(dom->div[div][0]))
4286 0 : return -1;
4287 0 : if (isl_seq_first_non_zero(dom->div[div] + 2 + d_dim, dom->n_div) != -1)
4288 0 : return -1;
4289 :
4290 0 : for (i = 0; i < bmap->n_div; ++i) {
4291 0 : if (isl_int_is_zero(bmap->div[i][0]))
4292 0 : continue;
4293 0 : if (isl_seq_first_non_zero(bmap->div[i] + 2 + d_dim,
4294 0 : (b_dim - d_dim) + bmap->n_div) != -1)
4295 0 : continue;
4296 0 : if (isl_seq_eq(bmap->div[i], dom->div[div], 2 + d_dim))
4297 0 : return i;
4298 : }
4299 0 : return -1;
4300 : }
4301 :
4302 : /* The correspondence between the variables in the main tableau,
4303 : * the context tableau, and the input map and domain is as follows.
4304 : * The first n_param and the last n_div variables of the main tableau
4305 : * form the variables of the context tableau.
4306 : * In the basic map, these n_param variables correspond to the
4307 : * parameters and the input dimensions. In the domain, they correspond
4308 : * to the parameters and the set dimensions.
4309 : * The n_div variables correspond to the integer divisions in the domain.
4310 : * To ensure that everything lines up, we may need to copy some of the
4311 : * integer divisions of the domain to the map. These have to be placed
4312 : * in the same order as those in the context and they have to be placed
4313 : * after any other integer divisions that the map may have.
4314 : * This function performs the required reordering.
4315 : */
4316 0 : static __isl_give isl_basic_map *align_context_divs(
4317 : __isl_take isl_basic_map *bmap, __isl_keep isl_basic_set *dom)
4318 : {
4319 : int i;
4320 0 : int common = 0;
4321 : int other;
4322 :
4323 0 : for (i = 0; i < dom->n_div; ++i)
4324 0 : if (find_context_div(bmap, dom, i) != -1)
4325 0 : common++;
4326 0 : other = bmap->n_div - common;
4327 0 : if (dom->n_div - common > 0) {
4328 0 : bmap = isl_basic_map_extend_space(bmap, isl_space_copy(bmap->dim),
4329 0 : dom->n_div - common, 0, 0);
4330 0 : if (!bmap)
4331 0 : return NULL;
4332 : }
4333 0 : for (i = 0; i < dom->n_div; ++i) {
4334 0 : int pos = find_context_div(bmap, dom, i);
4335 0 : if (pos < 0) {
4336 0 : pos = isl_basic_map_alloc_div(bmap);
4337 0 : if (pos < 0)
4338 0 : goto error;
4339 0 : isl_int_set_si(bmap->div[pos][0], 0);
4340 : }
4341 0 : if (pos != other + i)
4342 0 : isl_basic_map_swap_div(bmap, pos, other + i);
4343 : }
4344 0 : return bmap;
4345 : error:
4346 0 : isl_basic_map_free(bmap);
4347 0 : return NULL;
4348 : }
4349 :
4350 : /* Base case of isl_tab_basic_map_partial_lexopt, after removing
4351 : * some obvious symmetries.
4352 : *
4353 : * We make sure the divs in the domain are properly ordered,
4354 : * because they will be added one by one in the given order
4355 : * during the construction of the solution map.
4356 : * Furthermore, make sure that the known integer divisions
4357 : * appear before any unknown integer division because the solution
4358 : * may depend on the known integer divisions, while anything that
4359 : * depends on any variable starting from the first unknown integer
4360 : * division is ignored in sol_pma_add.
4361 : */
4362 0 : static struct isl_sol *basic_map_partial_lexopt_base_sol(
4363 : __isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom,
4364 : __isl_give isl_set **empty, int max,
4365 : struct isl_sol *(*init)(__isl_keep isl_basic_map *bmap,
4366 : __isl_take isl_basic_set *dom, int track_empty, int max))
4367 : {
4368 : struct isl_tab *tab;
4369 0 : struct isl_sol *sol = NULL;
4370 : struct isl_context *context;
4371 :
4372 0 : if (dom->n_div) {
4373 0 : dom = isl_basic_set_sort_divs(dom);
4374 0 : bmap = align_context_divs(bmap, dom);
4375 : }
4376 0 : sol = init(bmap, dom, !!empty, max);
4377 0 : if (!sol)
4378 0 : goto error;
4379 :
4380 0 : context = sol->context;
4381 0 : if (isl_basic_set_plain_is_empty(context->op->peek_basic_set(context)))
4382 : /* nothing */;
4383 0 : else if (isl_basic_map_plain_is_empty(bmap)) {
4384 0 : if (sol->add_empty)
4385 0 : sol->add_empty(sol,
4386 0 : isl_basic_set_copy(context->op->peek_basic_set(context)));
4387 : } else {
4388 0 : tab = tab_for_lexmin(bmap,
4389 0 : context->op->peek_basic_set(context), 1, max);
4390 0 : tab = context->op->detect_nonnegative_parameters(context, tab);
4391 0 : find_solutions_main(sol, tab);
4392 : }
4393 0 : if (sol->error)
4394 0 : goto error;
4395 :
4396 0 : isl_basic_map_free(bmap);
4397 0 : return sol;
4398 : error:
4399 0 : sol_free(sol);
4400 0 : isl_basic_map_free(bmap);
4401 0 : return NULL;
4402 : }
4403 :
4404 : /* Base case of isl_tab_basic_map_partial_lexopt, after removing
4405 : * some obvious symmetries.
4406 : *
4407 : * We call basic_map_partial_lexopt_base_sol and extract the results.
4408 : */
4409 0 : static __isl_give isl_map *basic_map_partial_lexopt_base(
4410 : __isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom,
4411 : __isl_give isl_set **empty, int max)
4412 : {
4413 0 : isl_map *result = NULL;
4414 : struct isl_sol *sol;
4415 : struct isl_sol_map *sol_map;
4416 :
4417 0 : sol = basic_map_partial_lexopt_base_sol(bmap, dom, empty, max,
4418 : &sol_map_init);
4419 0 : if (!sol)
4420 0 : return NULL;
4421 0 : sol_map = (struct isl_sol_map *) sol;
4422 :
4423 0 : result = isl_map_copy(sol_map->map);
4424 0 : if (empty)
4425 0 : *empty = isl_set_copy(sol_map->empty);
4426 0 : sol_free(&sol_map->sol);
4427 0 : return result;
4428 : }
4429 :
4430 : /* Return a count of the number of occurrences of the "n" first
4431 : * variables in the inequality constraints of "bmap".
4432 : */
4433 0 : static __isl_give int *count_occurrences(__isl_keep isl_basic_map *bmap,
4434 : int n)
4435 : {
4436 : int i, j;
4437 : isl_ctx *ctx;
4438 : int *occurrences;
4439 :
4440 0 : if (!bmap)
4441 0 : return NULL;
4442 0 : ctx = isl_basic_map_get_ctx(bmap);
4443 0 : occurrences = isl_calloc_array(ctx, int, n);
4444 0 : if (!occurrences)
4445 0 : return NULL;
4446 :
4447 0 : for (i = 0; i < bmap->n_ineq; ++i) {
4448 0 : for (j = 0; j < n; ++j) {
4449 0 : if (!isl_int_is_zero(bmap->ineq[i][1 + j]))
4450 0 : occurrences[j]++;
4451 : }
4452 : }
4453 :
4454 0 : return occurrences;
4455 : }
4456 :
4457 : /* Do all of the "n" variables with non-zero coefficients in "c"
4458 : * occur in exactly a single constraint.
4459 : * "occurrences" is an array of length "n" containing the number
4460 : * of occurrences of each of the variables in the inequality constraints.
4461 : */
4462 0 : static int single_occurrence(int n, isl_int *c, int *occurrences)
4463 : {
4464 : int i;
4465 :
4466 0 : for (i = 0; i < n; ++i) {
4467 0 : if (isl_int_is_zero(c[i]))
4468 0 : continue;
4469 0 : if (occurrences[i] != 1)
4470 0 : return 0;
4471 : }
4472 :
4473 0 : return 1;
4474 : }
4475 :
4476 : /* Do all of the "n" initial variables that occur in inequality constraint
4477 : * "ineq" of "bmap" only occur in that constraint?
4478 : */
4479 0 : static int all_single_occurrence(__isl_keep isl_basic_map *bmap, int ineq,
4480 : int n)
4481 : {
4482 : int i, j;
4483 :
4484 0 : for (i = 0; i < n; ++i) {
4485 0 : if (isl_int_is_zero(bmap->ineq[ineq][1 + i]))
4486 0 : continue;
4487 0 : for (j = 0; j < bmap->n_ineq; ++j) {
4488 0 : if (j == ineq)
4489 0 : continue;
4490 0 : if (!isl_int_is_zero(bmap->ineq[j][1 + i]))
4491 0 : return 0;
4492 : }
4493 : }
4494 :
4495 0 : return 1;
4496 : }
4497 :
4498 : /* Structure used during detection of parallel constraints.
4499 : * n_in: number of "input" variables: isl_dim_param + isl_dim_in
4500 : * n_out: number of "output" variables: isl_dim_out + isl_dim_div
4501 : * val: the coefficients of the output variables
4502 : */
4503 : struct isl_constraint_equal_info {
4504 : unsigned n_in;
4505 : unsigned n_out;
4506 : isl_int *val;
4507 : };
4508 :
4509 : /* Check whether the coefficients of the output variables
4510 : * of the constraint in "entry" are equal to info->val.
4511 : */
4512 0 : static int constraint_equal(const void *entry, const void *val)
4513 : {
4514 0 : isl_int **row = (isl_int **)entry;
4515 0 : const struct isl_constraint_equal_info *info = val;
4516 :
4517 0 : return isl_seq_eq((*row) + 1 + info->n_in, info->val, info->n_out);
4518 : }
4519 :
4520 : /* Check whether "bmap" has a pair of constraints that have
4521 : * the same coefficients for the output variables.
4522 : * Note that the coefficients of the existentially quantified
4523 : * variables need to be zero since the existentially quantified
4524 : * of the result are usually not the same as those of the input.
4525 : * Furthermore, check that each of the input variables that occur
4526 : * in those constraints does not occur in any other constraint.
4527 : * If so, return true and return the row indices of the two constraints
4528 : * in *first and *second.
4529 : */
4530 0 : static isl_bool parallel_constraints(__isl_keep isl_basic_map *bmap,
4531 : int *first, int *second)
4532 : {
4533 : int i;
4534 : isl_ctx *ctx;
4535 0 : int *occurrences = NULL;
4536 0 : struct isl_hash_table *table = NULL;
4537 : struct isl_hash_table_entry *entry;
4538 : struct isl_constraint_equal_info info;
4539 : unsigned n_out;
4540 : unsigned n_div;
4541 :
4542 0 : ctx = isl_basic_map_get_ctx(bmap);
4543 0 : table = isl_hash_table_alloc(ctx, bmap->n_ineq);
4544 0 : if (!table)
4545 0 : goto error;
4546 :
4547 0 : info.n_in = isl_basic_map_dim(bmap, isl_dim_param) +
4548 0 : isl_basic_map_dim(bmap, isl_dim_in);
4549 0 : occurrences = count_occurrences(bmap, info.n_in);
4550 0 : if (info.n_in && !occurrences)
4551 0 : goto error;
4552 0 : n_out = isl_basic_map_dim(bmap, isl_dim_out);
4553 0 : n_div = isl_basic_map_dim(bmap, isl_dim_div);
4554 0 : info.n_out = n_out + n_div;
4555 0 : for (i = 0; i < bmap->n_ineq; ++i) {
4556 : uint32_t hash;
4557 :
4558 0 : info.val = bmap->ineq[i] + 1 + info.n_in;
4559 0 : if (isl_seq_first_non_zero(info.val, n_out) < 0)
4560 0 : continue;
4561 0 : if (isl_seq_first_non_zero(info.val + n_out, n_div) >= 0)
4562 0 : continue;
4563 0 : if (!single_occurrence(info.n_in, bmap->ineq[i] + 1,
4564 : occurrences))
4565 0 : continue;
4566 0 : hash = isl_seq_get_hash(info.val, info.n_out);
4567 0 : entry = isl_hash_table_find(ctx, table, hash,
4568 : constraint_equal, &info, 1);
4569 0 : if (!entry)
4570 0 : goto error;
4571 0 : if (entry->data)
4572 0 : break;
4573 0 : entry->data = &bmap->ineq[i];
4574 : }
4575 :
4576 0 : if (i < bmap->n_ineq) {
4577 0 : *first = ((isl_int **)entry->data) - bmap->ineq;
4578 0 : *second = i;
4579 : }
4580 :
4581 0 : isl_hash_table_free(ctx, table);
4582 0 : free(occurrences);
4583 :
4584 0 : return i < bmap->n_ineq;
4585 : error:
4586 0 : isl_hash_table_free(ctx, table);
4587 0 : free(occurrences);
4588 0 : return isl_bool_error;
4589 : }
4590 :
4591 : /* Given a set of upper bounds in "var", add constraints to "bset"
4592 : * that make the i-th bound smallest.
4593 : *
4594 : * In particular, if there are n bounds b_i, then add the constraints
4595 : *
4596 : * b_i <= b_j for j > i
4597 : * b_i < b_j for j < i
4598 : */
4599 0 : static __isl_give isl_basic_set *select_minimum(__isl_take isl_basic_set *bset,
4600 : __isl_keep isl_mat *var, int i)
4601 : {
4602 : isl_ctx *ctx;
4603 : int j, k;
4604 :
4605 0 : ctx = isl_mat_get_ctx(var);
4606 :
4607 0 : for (j = 0; j < var->n_row; ++j) {
4608 0 : if (j == i)
4609 0 : continue;
4610 0 : k = isl_basic_set_alloc_inequality(bset);
4611 0 : if (k < 0)
4612 0 : goto error;
4613 0 : isl_seq_combine(bset->ineq[k], ctx->one, var->row[j],
4614 0 : ctx->negone, var->row[i], var->n_col);
4615 0 : isl_int_set_si(bset->ineq[k][var->n_col], 0);
4616 0 : if (j < i)
4617 0 : isl_int_sub_ui(bset->ineq[k][0], bset->ineq[k][0], 1);
4618 : }
4619 :
4620 0 : bset = isl_basic_set_finalize(bset);
4621 :
4622 0 : return bset;
4623 : error:
4624 0 : isl_basic_set_free(bset);
4625 0 : return NULL;
4626 : }
4627 :
4628 : /* Given a set of upper bounds on the last "input" variable m,
4629 : * construct a set that assigns the minimal upper bound to m, i.e.,
4630 : * construct a set that divides the space into cells where one
4631 : * of the upper bounds is smaller than all the others and assign
4632 : * this upper bound to m.
4633 : *
4634 : * In particular, if there are n bounds b_i, then the result
4635 : * consists of n basic sets, each one of the form
4636 : *
4637 : * m = b_i
4638 : * b_i <= b_j for j > i
4639 : * b_i < b_j for j < i
4640 : */
4641 0 : static __isl_give isl_set *set_minimum(__isl_take isl_space *dim,
4642 : __isl_take isl_mat *var)
4643 : {
4644 : int i, k;
4645 0 : isl_basic_set *bset = NULL;
4646 0 : isl_set *set = NULL;
4647 :
4648 0 : if (!dim || !var)
4649 : goto error;
4650 :
4651 0 : set = isl_set_alloc_space(isl_space_copy(dim),
4652 0 : var->n_row, ISL_SET_DISJOINT);
4653 :
4654 0 : for (i = 0; i < var->n_row; ++i) {
4655 0 : bset = isl_basic_set_alloc_space(isl_space_copy(dim), 0,
4656 0 : 1, var->n_row - 1);
4657 0 : k = isl_basic_set_alloc_equality(bset);
4658 0 : if (k < 0)
4659 0 : goto error;
4660 0 : isl_seq_cpy(bset->eq[k], var->row[i], var->n_col);
4661 0 : isl_int_set_si(bset->eq[k][var->n_col], -1);
4662 0 : bset = select_minimum(bset, var, i);
4663 0 : set = isl_set_add_basic_set(set, bset);
4664 : }
4665 :
4666 0 : isl_space_free(dim);
4667 0 : isl_mat_free(var);
4668 0 : return set;
4669 : error:
4670 0 : isl_basic_set_free(bset);
4671 0 : isl_set_free(set);
4672 0 : isl_space_free(dim);
4673 0 : isl_mat_free(var);
4674 0 : return NULL;
4675 : }
4676 :
4677 : /* Given that the last input variable of "bmap" represents the minimum
4678 : * of the bounds in "cst", check whether we need to split the domain
4679 : * based on which bound attains the minimum.
4680 : *
4681 : * A split is needed when the minimum appears in an integer division
4682 : * or in an equality. Otherwise, it is only needed if it appears in
4683 : * an upper bound that is different from the upper bounds on which it
4684 : * is defined.
4685 : */
4686 0 : static isl_bool need_split_basic_map(__isl_keep isl_basic_map *bmap,
4687 : __isl_keep isl_mat *cst)
4688 : {
4689 : int i, j;
4690 : unsigned total;
4691 : unsigned pos;
4692 :
4693 0 : pos = cst->n_col - 1;
4694 0 : total = isl_basic_map_dim(bmap, isl_dim_all);
4695 :
4696 0 : for (i = 0; i < bmap->n_div; ++i)
4697 0 : if (!isl_int_is_zero(bmap->div[i][2 + pos]))
4698 0 : return isl_bool_true;
4699 :
4700 0 : for (i = 0; i < bmap->n_eq; ++i)
4701 0 : if (!isl_int_is_zero(bmap->eq[i][1 + pos]))
4702 0 : return isl_bool_true;
4703 :
4704 0 : for (i = 0; i < bmap->n_ineq; ++i) {
4705 0 : if (isl_int_is_nonneg(bmap->ineq[i][1 + pos]))
4706 0 : continue;
4707 0 : if (!isl_int_is_negone(bmap->ineq[i][1 + pos]))
4708 0 : return isl_bool_true;
4709 0 : if (isl_seq_first_non_zero(bmap->ineq[i] + 1 + pos + 1,
4710 0 : total - pos - 1) >= 0)
4711 0 : return isl_bool_true;
4712 :
4713 0 : for (j = 0; j < cst->n_row; ++j)
4714 0 : if (isl_seq_eq(bmap->ineq[i], cst->row[j], cst->n_col))
4715 0 : break;
4716 0 : if (j >= cst->n_row)
4717 0 : return isl_bool_true;
4718 : }
4719 :
4720 0 : return isl_bool_false;
4721 : }
4722 :
4723 : /* Given that the last set variable of "bset" represents the minimum
4724 : * of the bounds in "cst", check whether we need to split the domain
4725 : * based on which bound attains the minimum.
4726 : *
4727 : * We simply call need_split_basic_map here. This is safe because
4728 : * the position of the minimum is computed from "cst" and not
4729 : * from "bmap".
4730 : */
4731 0 : static isl_bool need_split_basic_set(__isl_keep isl_basic_set *bset,
4732 : __isl_keep isl_mat *cst)
4733 : {
4734 0 : return need_split_basic_map(bset_to_bmap(bset), cst);
4735 : }
4736 :
4737 : /* Given that the last set variable of "set" represents the minimum
4738 : * of the bounds in "cst", check whether we need to split the domain
4739 : * based on which bound attains the minimum.
4740 : */
4741 0 : static isl_bool need_split_set(__isl_keep isl_set *set, __isl_keep isl_mat *cst)
4742 : {
4743 : int i;
4744 :
4745 0 : for (i = 0; i < set->n; ++i) {
4746 : isl_bool split;
4747 :
4748 0 : split = need_split_basic_set(set->p[i], cst);
4749 0 : if (split < 0 || split)
4750 0 : return split;
4751 : }
4752 :
4753 0 : return isl_bool_false;
4754 : }
4755 :
4756 : /* Given a set of which the last set variable is the minimum
4757 : * of the bounds in "cst", split each basic set in the set
4758 : * in pieces where one of the bounds is (strictly) smaller than the others.
4759 : * This subdivision is given in "min_expr".
4760 : * The variable is subsequently projected out.
4761 : *
4762 : * We only do the split when it is needed.
4763 : * For example if the last input variable m = min(a,b) and the only
4764 : * constraints in the given basic set are lower bounds on m,
4765 : * i.e., l <= m = min(a,b), then we can simply project out m
4766 : * to obtain l <= a and l <= b, without having to split on whether
4767 : * m is equal to a or b.
4768 : */
4769 0 : static __isl_give isl_set *split(__isl_take isl_set *empty,
4770 : __isl_take isl_set *min_expr, __isl_take isl_mat *cst)
4771 : {
4772 : int n_in;
4773 : int i;
4774 : isl_space *dim;
4775 : isl_set *res;
4776 :
4777 0 : if (!empty || !min_expr || !cst)
4778 : goto error;
4779 :
4780 0 : n_in = isl_set_dim(empty, isl_dim_set);
4781 0 : dim = isl_set_get_space(empty);
4782 0 : dim = isl_space_drop_dims(dim, isl_dim_set, n_in - 1, 1);
4783 0 : res = isl_set_empty(dim);
4784 :
4785 0 : for (i = 0; i < empty->n; ++i) {
4786 : isl_bool split;
4787 : isl_set *set;
4788 :
4789 0 : set = isl_set_from_basic_set(isl_basic_set_copy(empty->p[i]));
4790 0 : split = need_split_basic_set(empty->p[i], cst);
4791 0 : if (split < 0)
4792 0 : set = isl_set_free(set);
4793 0 : else if (split)
4794 0 : set = isl_set_intersect(set, isl_set_copy(min_expr));
4795 0 : set = isl_set_remove_dims(set, isl_dim_set, n_in - 1, 1);
4796 :
4797 0 : res = isl_set_union_disjoint(res, set);
4798 : }
4799 :
4800 0 : isl_set_free(empty);
4801 0 : isl_set_free(min_expr);
4802 0 : isl_mat_free(cst);
4803 0 : return res;
4804 : error:
4805 0 : isl_set_free(empty);
4806 0 : isl_set_free(min_expr);
4807 0 : isl_mat_free(cst);
4808 0 : return NULL;
4809 : }
4810 :
4811 : /* Given a map of which the last input variable is the minimum
4812 : * of the bounds in "cst", split each basic set in the set
4813 : * in pieces where one of the bounds is (strictly) smaller than the others.
4814 : * This subdivision is given in "min_expr".
4815 : * The variable is subsequently projected out.
4816 : *
4817 : * The implementation is essentially the same as that of "split".
4818 : */
4819 0 : static __isl_give isl_map *split_domain(__isl_take isl_map *opt,
4820 : __isl_take isl_set *min_expr, __isl_take isl_mat *cst)
4821 : {
4822 : int n_in;
4823 : int i;
4824 : isl_space *dim;
4825 : isl_map *res;
4826 :
4827 0 : if (!opt || !min_expr || !cst)
4828 : goto error;
4829 :
4830 0 : n_in = isl_map_dim(opt, isl_dim_in);
4831 0 : dim = isl_map_get_space(opt);
4832 0 : dim = isl_space_drop_dims(dim, isl_dim_in, n_in - 1, 1);
4833 0 : res = isl_map_empty(dim);
4834 :
4835 0 : for (i = 0; i < opt->n; ++i) {
4836 : isl_map *map;
4837 : isl_bool split;
4838 :
4839 0 : map = isl_map_from_basic_map(isl_basic_map_copy(opt->p[i]));
4840 0 : split = need_split_basic_map(opt->p[i], cst);
4841 0 : if (split < 0)
4842 0 : map = isl_map_free(map);
4843 0 : else if (split)
4844 0 : map = isl_map_intersect_domain(map,
4845 : isl_set_copy(min_expr));
4846 0 : map = isl_map_remove_dims(map, isl_dim_in, n_in - 1, 1);
4847 :
4848 0 : res = isl_map_union_disjoint(res, map);
4849 : }
4850 :
4851 0 : isl_map_free(opt);
4852 0 : isl_set_free(min_expr);
4853 0 : isl_mat_free(cst);
4854 0 : return res;
4855 : error:
4856 0 : isl_map_free(opt);
4857 0 : isl_set_free(min_expr);
4858 0 : isl_mat_free(cst);
4859 0 : return NULL;
4860 : }
4861 :
4862 : static __isl_give isl_map *basic_map_partial_lexopt(
4863 : __isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom,
4864 : __isl_give isl_set **empty, int max);
4865 :
4866 : /* This function is called from basic_map_partial_lexopt_symm.
4867 : * The last variable of "bmap" and "dom" corresponds to the minimum
4868 : * of the bounds in "cst". "map_space" is the space of the original
4869 : * input relation (of basic_map_partial_lexopt_symm) and "set_space"
4870 : * is the space of the original domain.
4871 : *
4872 : * We recursively call basic_map_partial_lexopt and then plug in
4873 : * the definition of the minimum in the result.
4874 : */
4875 0 : static __isl_give isl_map *basic_map_partial_lexopt_symm_core(
4876 : __isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom,
4877 : __isl_give isl_set **empty, int max, __isl_take isl_mat *cst,
4878 : __isl_take isl_space *map_space, __isl_take isl_space *set_space)
4879 : {
4880 : isl_map *opt;
4881 : isl_set *min_expr;
4882 :
4883 0 : min_expr = set_minimum(isl_basic_set_get_space(dom), isl_mat_copy(cst));
4884 :
4885 0 : opt = basic_map_partial_lexopt(bmap, dom, empty, max);
4886 :
4887 0 : if (empty) {
4888 0 : *empty = split(*empty,
4889 : isl_set_copy(min_expr), isl_mat_copy(cst));
4890 0 : *empty = isl_set_reset_space(*empty, set_space);
4891 : }
4892 :
4893 0 : opt = split_domain(opt, min_expr, cst);
4894 0 : opt = isl_map_reset_space(opt, map_space);
4895 :
4896 0 : return opt;
4897 : }
4898 :
4899 : /* Extract a domain from "bmap" for the purpose of computing
4900 : * a lexicographic optimum.
4901 : *
4902 : * This function is only called when the caller wants to compute a full
4903 : * lexicographic optimum, i.e., without specifying a domain. In this case,
4904 : * the caller is not interested in the part of the domain space where
4905 : * there is no solution and the domain can be initialized to those constraints
4906 : * of "bmap" that only involve the parameters and the input dimensions.
4907 : * This relieves the parametric programming engine from detecting those
4908 : * inequalities and transferring them to the context. More importantly,
4909 : * it ensures that those inequalities are transferred first and not
4910 : * intermixed with inequalities that actually split the domain.
4911 : *
4912 : * If the caller does not require the absence of existentially quantified
4913 : * variables in the result (i.e., if ISL_OPT_QE is not set in "flags"),
4914 : * then the actual domain of "bmap" can be used. This ensures that
4915 : * the domain does not need to be split at all just to separate out
4916 : * pieces of the domain that do not have a solution from piece that do.
4917 : * This domain cannot be used in general because it may involve
4918 : * (unknown) existentially quantified variables which will then also
4919 : * appear in the solution.
4920 : */
4921 0 : static __isl_give isl_basic_set *extract_domain(__isl_keep isl_basic_map *bmap,
4922 : unsigned flags)
4923 : {
4924 : int n_div;
4925 : int n_out;
4926 :
4927 0 : n_div = isl_basic_map_dim(bmap, isl_dim_div);
4928 0 : n_out = isl_basic_map_dim(bmap, isl_dim_out);
4929 0 : bmap = isl_basic_map_copy(bmap);
4930 0 : if (ISL_FL_ISSET(flags, ISL_OPT_QE)) {
4931 0 : bmap = isl_basic_map_drop_constraints_involving_dims(bmap,
4932 : isl_dim_div, 0, n_div);
4933 0 : bmap = isl_basic_map_drop_constraints_involving_dims(bmap,
4934 : isl_dim_out, 0, n_out);
4935 : }
4936 0 : return isl_basic_map_domain(bmap);
4937 : }
4938 :
4939 : #undef TYPE
4940 : #define TYPE isl_map
4941 : #undef SUFFIX
4942 : #define SUFFIX
4943 : #include "isl_tab_lexopt_templ.c"
4944 :
4945 : /* Extract the subsequence of the sample value of "tab"
4946 : * starting at "pos" and of length "len".
4947 : */
4948 0 : static __isl_give isl_vec *extract_sample_sequence(struct isl_tab *tab,
4949 : int pos, int len)
4950 : {
4951 : int i;
4952 : isl_ctx *ctx;
4953 : isl_vec *v;
4954 :
4955 0 : ctx = isl_tab_get_ctx(tab);
4956 0 : v = isl_vec_alloc(ctx, len);
4957 0 : if (!v)
4958 0 : return NULL;
4959 0 : for (i = 0; i < len; ++i) {
4960 0 : if (!tab->var[pos + i].is_row) {
4961 0 : isl_int_set_si(v->el[i], 0);
4962 : } else {
4963 : int row;
4964 :
4965 0 : row = tab->var[pos + i].index;
4966 0 : isl_int_divexact(v->el[i], tab->mat->row[row][1],
4967 : tab->mat->row[row][0]);
4968 : }
4969 : }
4970 :
4971 0 : return v;
4972 : }
4973 :
4974 : /* Check if the sequence of variables starting at "pos"
4975 : * represents a trivial solution according to "trivial".
4976 : * That is, is the result of applying "trivial" to this sequence
4977 : * equal to the zero vector?
4978 : */
4979 0 : static isl_bool region_is_trivial(struct isl_tab *tab, int pos,
4980 : __isl_keep isl_mat *trivial)
4981 : {
4982 : int n, len;
4983 : isl_vec *v;
4984 : isl_bool is_trivial;
4985 :
4986 0 : if (!trivial)
4987 0 : return isl_bool_error;
4988 :
4989 0 : n = isl_mat_rows(trivial);
4990 0 : if (n == 0)
4991 0 : return isl_bool_false;
4992 :
4993 0 : len = isl_mat_cols(trivial);
4994 0 : v = extract_sample_sequence(tab, pos, len);
4995 0 : v = isl_mat_vec_product(isl_mat_copy(trivial), v);
4996 0 : is_trivial = isl_vec_is_zero(v);
4997 0 : isl_vec_free(v);
4998 :
4999 0 : return is_trivial;
5000 : }
5001 :
5002 : /* Global internal data for isl_tab_basic_set_non_trivial_lexmin.
5003 : *
5004 : * "n_op" is the number of initial coordinates to optimize,
5005 : * as passed to isl_tab_basic_set_non_trivial_lexmin.
5006 : * "region" is the "n_region"-sized array of regions passed
5007 : * to isl_tab_basic_set_non_trivial_lexmin.
5008 : *
5009 : * "tab" is the tableau that corresponds to the ILP problem.
5010 : * "local" is an array of local data structure, one for each
5011 : * (potential) level of the backtracking procedure of
5012 : * isl_tab_basic_set_non_trivial_lexmin.
5013 : * "v" is a pre-allocated vector that can be used for adding
5014 : * constraints to the tableau.
5015 : *
5016 : * "sol" contains the best solution found so far.
5017 : * It is initialized to a vector of size zero.
5018 : */
5019 : struct isl_lexmin_data {
5020 : int n_op;
5021 : int n_region;
5022 : struct isl_trivial_region *region;
5023 :
5024 : struct isl_tab *tab;
5025 : struct isl_local_region *local;
5026 : isl_vec *v;
5027 :
5028 : isl_vec *sol;
5029 : };
5030 :
5031 : /* Return the index of the first trivial region, "n_region" if all regions
5032 : * are non-trivial or -1 in case of error.
5033 : */
5034 0 : static int first_trivial_region(struct isl_lexmin_data *data)
5035 : {
5036 : int i;
5037 :
5038 0 : for (i = 0; i < data->n_region; ++i) {
5039 : isl_bool trivial;
5040 0 : trivial = region_is_trivial(data->tab, data->region[i].pos,
5041 0 : data->region[i].trivial);
5042 0 : if (trivial < 0)
5043 0 : return -1;
5044 0 : if (trivial)
5045 0 : return i;
5046 : }
5047 :
5048 0 : return data->n_region;
5049 : }
5050 :
5051 : /* Check if the solution is optimal, i.e., whether the first
5052 : * n_op entries are zero.
5053 : */
5054 0 : static int is_optimal(__isl_keep isl_vec *sol, int n_op)
5055 : {
5056 : int i;
5057 :
5058 0 : for (i = 0; i < n_op; ++i)
5059 0 : if (!isl_int_is_zero(sol->el[1 + i]))
5060 0 : return 0;
5061 0 : return 1;
5062 : }
5063 :
5064 : /* Add constraints to "tab" that ensure that any solution is significantly
5065 : * better than that represented by "sol". That is, find the first
5066 : * relevant (within first n_op) non-zero coefficient and force it (along
5067 : * with all previous coefficients) to be zero.
5068 : * If the solution is already optimal (all relevant coefficients are zero),
5069 : * then just mark the table as empty.
5070 : * "n_zero" is the number of coefficients that have been forced zero
5071 : * by previous calls to this function at the same level.
5072 : * Return the updated number of forced zero coefficients or -1 on error.
5073 : *
5074 : * This function assumes that at least 2 * (n_op - n_zero) more rows and
5075 : * at least 2 * (n_op - n_zero) more elements in the constraint array
5076 : * are available in the tableau.
5077 : */
5078 0 : static int force_better_solution(struct isl_tab *tab,
5079 : __isl_keep isl_vec *sol, int n_op, int n_zero)
5080 : {
5081 : int i, n;
5082 : isl_ctx *ctx;
5083 0 : isl_vec *v = NULL;
5084 :
5085 0 : if (!sol)
5086 0 : return -1;
5087 :
5088 0 : for (i = n_zero; i < n_op; ++i)
5089 0 : if (!isl_int_is_zero(sol->el[1 + i]))
5090 : break;
5091 :
5092 0 : if (i == n_op) {
5093 0 : if (isl_tab_mark_empty(tab) < 0)
5094 0 : return -1;
5095 0 : return n_op;
5096 : }
5097 :
5098 0 : ctx = isl_vec_get_ctx(sol);
5099 0 : v = isl_vec_alloc(ctx, 1 + tab->n_var);
5100 0 : if (!v)
5101 0 : return -1;
5102 :
5103 0 : n = i + 1;
5104 0 : for (; i >= n_zero; --i) {
5105 0 : v = isl_vec_clr(v);
5106 0 : isl_int_set_si(v->el[1 + i], -1);
5107 0 : if (add_lexmin_eq(tab, v->el) < 0)
5108 0 : goto error;
5109 : }
5110 :
5111 0 : isl_vec_free(v);
5112 0 : return n;
5113 : error:
5114 0 : isl_vec_free(v);
5115 0 : return -1;
5116 : }
5117 :
5118 : /* Fix triviality direction "dir" of the given region to zero.
5119 : *
5120 : * This function assumes that at least two more rows and at least
5121 : * two more elements in the constraint array are available in the tableau.
5122 : */
5123 0 : static isl_stat fix_zero(struct isl_tab *tab, struct isl_trivial_region *region,
5124 : int dir, struct isl_lexmin_data *data)
5125 : {
5126 : int len;
5127 :
5128 0 : data->v = isl_vec_clr(data->v);
5129 0 : if (!data->v)
5130 0 : return isl_stat_error;
5131 0 : len = isl_mat_cols(region->trivial);
5132 0 : isl_seq_cpy(data->v->el + 1 + region->pos, region->trivial->row[dir],
5133 : len);
5134 0 : if (add_lexmin_eq(tab, data->v->el) < 0)
5135 0 : return isl_stat_error;
5136 :
5137 0 : return isl_stat_ok;
5138 : }
5139 :
5140 : /* This function selects case "side" for non-triviality region "region",
5141 : * assuming all the equality constraints have been imposed already.
5142 : * In particular, the triviality direction side/2 is made positive
5143 : * if side is even and made negative if side is odd.
5144 : *
5145 : * This function assumes that at least one more row and at least
5146 : * one more element in the constraint array are available in the tableau.
5147 : */
5148 0 : static struct isl_tab *pos_neg(struct isl_tab *tab,
5149 : struct isl_trivial_region *region,
5150 : int side, struct isl_lexmin_data *data)
5151 : {
5152 : int len;
5153 :
5154 0 : data->v = isl_vec_clr(data->v);
5155 0 : if (!data->v)
5156 0 : goto error;
5157 0 : isl_int_set_si(data->v->el[0], -1);
5158 0 : len = isl_mat_cols(region->trivial);
5159 0 : if (side % 2 == 0)
5160 0 : isl_seq_cpy(data->v->el + 1 + region->pos,
5161 0 : region->trivial->row[side / 2], len);
5162 : else
5163 0 : isl_seq_neg(data->v->el + 1 + region->pos,
5164 0 : region->trivial->row[side / 2], len);
5165 0 : return add_lexmin_ineq(tab, data->v->el);
5166 : error:
5167 0 : isl_tab_free(tab);
5168 0 : return NULL;
5169 : }
5170 :
5171 : /* Local data at each level of the backtracking procedure of
5172 : * isl_tab_basic_set_non_trivial_lexmin.
5173 : *
5174 : * "update" is set if a solution has been found in the current case
5175 : * of this level, such that a better solution needs to be enforced
5176 : * in the next case.
5177 : * "n_zero" is the number of initial coordinates that have already
5178 : * been forced to be zero at this level.
5179 : * "region" is the non-triviality region considered at this level.
5180 : * "side" is the index of the current case at this level.
5181 : * "n" is the number of triviality directions.
5182 : * "snap" is a snapshot of the tableau holding a state that needs
5183 : * to be satisfied by all subsequent cases.
5184 : */
5185 : struct isl_local_region {
5186 : int update;
5187 : int n_zero;
5188 : int region;
5189 : int side;
5190 : int n;
5191 : struct isl_tab_undo *snap;
5192 : };
5193 :
5194 : /* Initialize the global data structure "data" used while solving
5195 : * the ILP problem "bset".
5196 : */
5197 0 : static isl_stat init_lexmin_data(struct isl_lexmin_data *data,
5198 : __isl_keep isl_basic_set *bset)
5199 : {
5200 : isl_ctx *ctx;
5201 :
5202 0 : ctx = isl_basic_set_get_ctx(bset);
5203 :
5204 0 : data->tab = tab_for_lexmin(bset, NULL, 0, 0);
5205 0 : if (!data->tab)
5206 0 : return isl_stat_error;
5207 :
5208 0 : data->v = isl_vec_alloc(ctx, 1 + data->tab->n_var);
5209 0 : if (!data->v)
5210 0 : return isl_stat_error;
5211 0 : data->local = isl_calloc_array(ctx, struct isl_local_region,
5212 : data->n_region);
5213 0 : if (data->n_region && !data->local)
5214 0 : return isl_stat_error;
5215 :
5216 0 : data->sol = isl_vec_alloc(ctx, 0);
5217 :
5218 0 : return isl_stat_ok;
5219 : }
5220 :
5221 : /* Mark all outer levels as requiring a better solution
5222 : * in the next cases.
5223 : */
5224 0 : static void update_outer_levels(struct isl_lexmin_data *data, int level)
5225 : {
5226 : int i;
5227 :
5228 0 : for (i = 0; i < level; ++i)
5229 0 : data->local[i].update = 1;
5230 0 : }
5231 :
5232 : /* Initialize "local" to refer to region "region" and
5233 : * to initiate processing at this level.
5234 : */
5235 0 : static void init_local_region(struct isl_local_region *local, int region,
5236 : struct isl_lexmin_data *data)
5237 : {
5238 0 : local->n = isl_mat_rows(data->region[region].trivial);
5239 0 : local->region = region;
5240 0 : local->side = 0;
5241 0 : local->update = 0;
5242 0 : local->n_zero = 0;
5243 0 : }
5244 :
5245 : /* What to do next after entering a level of the backtracking procedure.
5246 : *
5247 : * error: some error has occurred; abort
5248 : * done: an optimal solution has been found; stop search
5249 : * backtrack: backtrack to the previous level
5250 : * handle: add the constraints for the current level and
5251 : * move to the next level
5252 : */
5253 : enum isl_next {
5254 : isl_next_error = -1,
5255 : isl_next_done,
5256 : isl_next_backtrack,
5257 : isl_next_handle,
5258 : };
5259 :
5260 : /* Have all cases of the current region been considered?
5261 : * If there are n directions, then there are 2n cases.
5262 : *
5263 : * The constraints in the current tableau are imposed
5264 : * in all subsequent cases. This means that if the current
5265 : * tableau is empty, then none of those cases should be considered
5266 : * anymore and all cases have effectively been considered.
5267 : */
5268 0 : static int finished_all_cases(struct isl_local_region *local,
5269 : struct isl_lexmin_data *data)
5270 : {
5271 0 : if (data->tab->empty)
5272 0 : return 1;
5273 0 : return local->side >= 2 * local->n;
5274 : }
5275 :
5276 : /* Enter level "level" of the backtracking search and figure out
5277 : * what to do next. "init" is set if the level was entered
5278 : * from a higher level and needs to be initialized.
5279 : * Otherwise, the level is entered as a result of backtracking and
5280 : * the tableau needs to be restored to a position that can
5281 : * be used for the next case at this level.
5282 : * The snapshot is assumed to have been saved in the previous case,
5283 : * before the constraints specific to that case were added.
5284 : *
5285 : * In the initialization case, the local region is initialized
5286 : * to point to the first violated region.
5287 : * If the constraints of all regions are satisfied by the current
5288 : * sample of the tableau, then tell the caller to continue looking
5289 : * for a better solution or to stop searching if an optimal solution
5290 : * has been found.
5291 : *
5292 : * If the tableau is empty or if all cases at the current level
5293 : * have been considered, then the caller needs to backtrack as well.
5294 : */
5295 0 : static enum isl_next enter_level(int level, int init,
5296 : struct isl_lexmin_data *data)
5297 : {
5298 0 : struct isl_local_region *local = &data->local[level];
5299 :
5300 0 : if (init) {
5301 : int r;
5302 :
5303 0 : data->tab = cut_to_integer_lexmin(data->tab, CUT_ONE);
5304 0 : if (!data->tab)
5305 0 : return isl_next_error;
5306 0 : if (data->tab->empty)
5307 0 : return isl_next_backtrack;
5308 0 : r = first_trivial_region(data);
5309 0 : if (r < 0)
5310 0 : return isl_next_error;
5311 0 : if (r == data->n_region) {
5312 0 : update_outer_levels(data, level);
5313 0 : isl_vec_free(data->sol);
5314 0 : data->sol = isl_tab_get_sample_value(data->tab);
5315 0 : if (!data->sol)
5316 0 : return isl_next_error;
5317 0 : if (is_optimal(data->sol, data->n_op))
5318 0 : return isl_next_done;
5319 0 : return isl_next_backtrack;
5320 : }
5321 0 : if (level >= data->n_region)
5322 0 : isl_die(isl_vec_get_ctx(data->v), isl_error_internal,
5323 : "nesting level too deep",
5324 : return isl_next_error);
5325 0 : init_local_region(local, r, data);
5326 0 : if (isl_tab_extend_cons(data->tab,
5327 0 : 2 * local->n + 2 * data->n_op) < 0)
5328 0 : return isl_next_error;
5329 : } else {
5330 0 : if (isl_tab_rollback(data->tab, local->snap) < 0)
5331 0 : return isl_next_error;
5332 : }
5333 :
5334 0 : if (finished_all_cases(local, data))
5335 0 : return isl_next_backtrack;
5336 0 : return isl_next_handle;
5337 : }
5338 :
5339 : /* If a solution has been found in the previous case at this level
5340 : * (marked by local->update being set), then add constraints
5341 : * that enforce a better solution in the present and all following cases.
5342 : * The constraints only need to be imposed once because they are
5343 : * included in the snapshot (taken in pick_side) that will be used in
5344 : * subsequent cases.
5345 : */
5346 0 : static isl_stat better_next_side(struct isl_local_region *local,
5347 : struct isl_lexmin_data *data)
5348 : {
5349 0 : if (!local->update)
5350 0 : return isl_stat_ok;
5351 :
5352 0 : local->n_zero = force_better_solution(data->tab,
5353 : data->sol, data->n_op, local->n_zero);
5354 0 : if (local->n_zero < 0)
5355 0 : return isl_stat_error;
5356 :
5357 0 : local->update = 0;
5358 :
5359 0 : return isl_stat_ok;
5360 : }
5361 :
5362 : /* Add constraints to data->tab that select the current case (local->side)
5363 : * at the current level.
5364 : *
5365 : * If the linear combinations v should not be zero, then the cases are
5366 : * v_0 >= 1
5367 : * v_0 <= -1
5368 : * v_0 = 0 and v_1 >= 1
5369 : * v_0 = 0 and v_1 <= -1
5370 : * v_0 = 0 and v_1 = 0 and v_2 >= 1
5371 : * v_0 = 0 and v_1 = 0 and v_2 <= -1
5372 : * ...
5373 : * in this order.
5374 : *
5375 : * A snapshot is taken after the equality constraint (if any) has been added
5376 : * such that the next case can start off from this position.
5377 : * The rollback to this position is performed in enter_level.
5378 : */
5379 0 : static isl_stat pick_side(struct isl_local_region *local,
5380 : struct isl_lexmin_data *data)
5381 : {
5382 : struct isl_trivial_region *region;
5383 : int side, base;
5384 :
5385 0 : region = &data->region[local->region];
5386 0 : side = local->side;
5387 0 : base = 2 * (side/2);
5388 :
5389 0 : if (side == base && base >= 2 &&
5390 0 : fix_zero(data->tab, region, base / 2 - 1, data) < 0)
5391 0 : return isl_stat_error;
5392 :
5393 0 : local->snap = isl_tab_snap(data->tab);
5394 0 : if (isl_tab_push_basis(data->tab) < 0)
5395 0 : return isl_stat_error;
5396 :
5397 0 : data->tab = pos_neg(data->tab, region, side, data);
5398 0 : if (!data->tab)
5399 0 : return isl_stat_error;
5400 0 : return isl_stat_ok;
5401 : }
5402 :
5403 : /* Free the memory associated to "data".
5404 : */
5405 0 : static void clear_lexmin_data(struct isl_lexmin_data *data)
5406 : {
5407 0 : free(data->local);
5408 0 : isl_vec_free(data->v);
5409 0 : isl_tab_free(data->tab);
5410 0 : }
5411 :
5412 : /* Return the lexicographically smallest non-trivial solution of the
5413 : * given ILP problem.
5414 : *
5415 : * All variables are assumed to be non-negative.
5416 : *
5417 : * n_op is the number of initial coordinates to optimize.
5418 : * That is, once a solution has been found, we will only continue looking
5419 : * for solutions that result in significantly better values for those
5420 : * initial coordinates. That is, we only continue looking for solutions
5421 : * that increase the number of initial zeros in this sequence.
5422 : *
5423 : * A solution is non-trivial, if it is non-trivial on each of the
5424 : * specified regions. Each region represents a sequence of
5425 : * triviality directions on a sequence of variables that starts
5426 : * at a given position. A solution is non-trivial on such a region if
5427 : * at least one of the triviality directions is non-zero
5428 : * on that sequence of variables.
5429 : *
5430 : * Whenever a conflict is encountered, all constraints involved are
5431 : * reported to the caller through a call to "conflict".
5432 : *
5433 : * We perform a simple branch-and-bound backtracking search.
5434 : * Each level in the search represents an initially trivial region
5435 : * that is forced to be non-trivial.
5436 : * At each level we consider 2 * n cases, where n
5437 : * is the number of triviality directions.
5438 : * In terms of those n directions v_i, we consider the cases
5439 : * v_0 >= 1
5440 : * v_0 <= -1
5441 : * v_0 = 0 and v_1 >= 1
5442 : * v_0 = 0 and v_1 <= -1
5443 : * v_0 = 0 and v_1 = 0 and v_2 >= 1
5444 : * v_0 = 0 and v_1 = 0 and v_2 <= -1
5445 : * ...
5446 : * in this order.
5447 : */
5448 0 : __isl_give isl_vec *isl_tab_basic_set_non_trivial_lexmin(
5449 : __isl_take isl_basic_set *bset, int n_op, int n_region,
5450 : struct isl_trivial_region *region,
5451 : int (*conflict)(int con, void *user), void *user)
5452 : {
5453 0 : struct isl_lexmin_data data = { n_op, n_region, region };
5454 : int level, init;
5455 :
5456 0 : if (!bset)
5457 0 : return NULL;
5458 :
5459 0 : if (init_lexmin_data(&data, bset) < 0)
5460 0 : goto error;
5461 0 : data.tab->conflict = conflict;
5462 0 : data.tab->conflict_user = user;
5463 :
5464 0 : level = 0;
5465 0 : init = 1;
5466 :
5467 0 : while (level >= 0) {
5468 : enum isl_next next;
5469 0 : struct isl_local_region *local = &data.local[level];
5470 :
5471 0 : next = enter_level(level, init, &data);
5472 0 : if (next < 0)
5473 0 : goto error;
5474 0 : if (next == isl_next_done)
5475 0 : break;
5476 0 : if (next == isl_next_backtrack) {
5477 0 : level--;
5478 0 : init = 0;
5479 0 : continue;
5480 : }
5481 :
5482 0 : if (better_next_side(local, &data) < 0)
5483 0 : goto error;
5484 0 : if (pick_side(local, &data) < 0)
5485 0 : goto error;
5486 :
5487 0 : local->side++;
5488 0 : level++;
5489 0 : init = 1;
5490 : }
5491 :
5492 0 : clear_lexmin_data(&data);
5493 0 : isl_basic_set_free(bset);
5494 :
5495 0 : return data.sol;
5496 : error:
5497 0 : clear_lexmin_data(&data);
5498 0 : isl_basic_set_free(bset);
5499 0 : isl_vec_free(data.sol);
5500 0 : return NULL;
5501 : }
5502 :
5503 : /* Wrapper for a tableau that is used for computing
5504 : * the lexicographically smallest rational point of a non-negative set.
5505 : * This point is represented by the sample value of "tab",
5506 : * unless "tab" is empty.
5507 : */
5508 : struct isl_tab_lexmin {
5509 : isl_ctx *ctx;
5510 : struct isl_tab *tab;
5511 : };
5512 :
5513 : /* Free "tl" and return NULL.
5514 : */
5515 0 : __isl_null isl_tab_lexmin *isl_tab_lexmin_free(__isl_take isl_tab_lexmin *tl)
5516 : {
5517 0 : if (!tl)
5518 0 : return NULL;
5519 0 : isl_ctx_deref(tl->ctx);
5520 0 : isl_tab_free(tl->tab);
5521 0 : free(tl);
5522 :
5523 0 : return NULL;
5524 : }
5525 :
5526 : /* Construct an isl_tab_lexmin for computing
5527 : * the lexicographically smallest rational point in "bset",
5528 : * assuming that all variables are non-negative.
5529 : */
5530 0 : __isl_give isl_tab_lexmin *isl_tab_lexmin_from_basic_set(
5531 : __isl_take isl_basic_set *bset)
5532 : {
5533 : isl_ctx *ctx;
5534 : isl_tab_lexmin *tl;
5535 :
5536 0 : if (!bset)
5537 0 : return NULL;
5538 :
5539 0 : ctx = isl_basic_set_get_ctx(bset);
5540 0 : tl = isl_calloc_type(ctx, struct isl_tab_lexmin);
5541 0 : if (!tl)
5542 0 : goto error;
5543 0 : tl->ctx = ctx;
5544 0 : isl_ctx_ref(ctx);
5545 0 : tl->tab = tab_for_lexmin(bset, NULL, 0, 0);
5546 0 : isl_basic_set_free(bset);
5547 0 : if (!tl->tab)
5548 0 : return isl_tab_lexmin_free(tl);
5549 0 : return tl;
5550 : error:
5551 0 : isl_basic_set_free(bset);
5552 0 : isl_tab_lexmin_free(tl);
5553 0 : return NULL;
5554 : }
5555 :
5556 : /* Return the dimension of the set represented by "tl".
5557 : */
5558 0 : int isl_tab_lexmin_dim(__isl_keep isl_tab_lexmin *tl)
5559 : {
5560 0 : return tl ? tl->tab->n_var : -1;
5561 : }
5562 :
5563 : /* Add the equality with coefficients "eq" to "tl", updating the optimal
5564 : * solution if needed.
5565 : * The equality is added as two opposite inequality constraints.
5566 : */
5567 0 : __isl_give isl_tab_lexmin *isl_tab_lexmin_add_eq(__isl_take isl_tab_lexmin *tl,
5568 : isl_int *eq)
5569 : {
5570 : unsigned n_var;
5571 :
5572 0 : if (!tl || !eq)
5573 0 : return isl_tab_lexmin_free(tl);
5574 :
5575 0 : if (isl_tab_extend_cons(tl->tab, 2) < 0)
5576 0 : return isl_tab_lexmin_free(tl);
5577 0 : n_var = tl->tab->n_var;
5578 0 : isl_seq_neg(eq, eq, 1 + n_var);
5579 0 : tl->tab = add_lexmin_ineq(tl->tab, eq);
5580 0 : isl_seq_neg(eq, eq, 1 + n_var);
5581 0 : tl->tab = add_lexmin_ineq(tl->tab, eq);
5582 :
5583 0 : if (!tl->tab)
5584 0 : return isl_tab_lexmin_free(tl);
5585 :
5586 0 : return tl;
5587 : }
5588 :
5589 : /* Add cuts to "tl" until the sample value reaches an integer value or
5590 : * until the result becomes empty.
5591 : */
5592 0 : __isl_give isl_tab_lexmin *isl_tab_lexmin_cut_to_integer(
5593 : __isl_take isl_tab_lexmin *tl)
5594 : {
5595 0 : if (!tl)
5596 0 : return NULL;
5597 0 : tl->tab = cut_to_integer_lexmin(tl->tab, CUT_ONE);
5598 0 : if (!tl->tab)
5599 0 : return isl_tab_lexmin_free(tl);
5600 0 : return tl;
5601 : }
5602 :
5603 : /* Return the lexicographically smallest rational point in the basic set
5604 : * from which "tl" was constructed.
5605 : * If the original input was empty, then return a zero-length vector.
5606 : */
5607 0 : __isl_give isl_vec *isl_tab_lexmin_get_solution(__isl_keep isl_tab_lexmin *tl)
5608 : {
5609 0 : if (!tl)
5610 0 : return NULL;
5611 0 : if (tl->tab->empty)
5612 0 : return isl_vec_alloc(tl->ctx, 0);
5613 : else
5614 0 : return isl_tab_get_sample_value(tl->tab);
5615 : }
5616 :
5617 : struct isl_sol_pma {
5618 : struct isl_sol sol;
5619 : isl_pw_multi_aff *pma;
5620 : isl_set *empty;
5621 : };
5622 :
5623 0 : static void sol_pma_free(struct isl_sol *sol)
5624 : {
5625 0 : struct isl_sol_pma *sol_pma = (struct isl_sol_pma *) sol;
5626 0 : isl_pw_multi_aff_free(sol_pma->pma);
5627 0 : isl_set_free(sol_pma->empty);
5628 0 : }
5629 :
5630 : /* This function is called for parts of the context where there is
5631 : * no solution, with "bset" corresponding to the context tableau.
5632 : * Simply add the basic set to the set "empty".
5633 : */
5634 0 : static void sol_pma_add_empty(struct isl_sol_pma *sol,
5635 : __isl_take isl_basic_set *bset)
5636 : {
5637 0 : if (!bset || !sol->empty)
5638 : goto error;
5639 :
5640 0 : sol->empty = isl_set_grow(sol->empty, 1);
5641 0 : bset = isl_basic_set_simplify(bset);
5642 0 : bset = isl_basic_set_finalize(bset);
5643 0 : sol->empty = isl_set_add_basic_set(sol->empty, bset);
5644 0 : if (!sol->empty)
5645 0 : sol->sol.error = 1;
5646 0 : return;
5647 : error:
5648 0 : isl_basic_set_free(bset);
5649 0 : sol->sol.error = 1;
5650 : }
5651 :
5652 : /* Given a basic set "dom" that represents the context and a tuple of
5653 : * affine expressions "maff" defined over this domain, construct
5654 : * an isl_pw_multi_aff with a single cell corresponding to "dom" and
5655 : * the affine expressions in "maff".
5656 : */
5657 0 : static void sol_pma_add(struct isl_sol_pma *sol,
5658 : __isl_take isl_basic_set *dom, __isl_take isl_multi_aff *maff)
5659 : {
5660 : isl_pw_multi_aff *pma;
5661 :
5662 0 : dom = isl_basic_set_simplify(dom);
5663 0 : dom = isl_basic_set_finalize(dom);
5664 0 : pma = isl_pw_multi_aff_alloc(isl_set_from_basic_set(dom), maff);
5665 0 : sol->pma = isl_pw_multi_aff_add_disjoint(sol->pma, pma);
5666 0 : if (!sol->pma)
5667 0 : sol->sol.error = 1;
5668 0 : }
5669 :
5670 0 : static void sol_pma_add_empty_wrap(struct isl_sol *sol,
5671 : __isl_take isl_basic_set *bset)
5672 : {
5673 0 : sol_pma_add_empty((struct isl_sol_pma *)sol, bset);
5674 0 : }
5675 :
5676 0 : static void sol_pma_add_wrap(struct isl_sol *sol,
5677 : __isl_take isl_basic_set *dom, __isl_take isl_multi_aff *ma)
5678 : {
5679 0 : sol_pma_add((struct isl_sol_pma *)sol, dom, ma);
5680 0 : }
5681 :
5682 : /* Construct an isl_sol_pma structure for accumulating the solution.
5683 : * If track_empty is set, then we also keep track of the parts
5684 : * of the context where there is no solution.
5685 : * If max is set, then we are solving a maximization, rather than
5686 : * a minimization problem, which means that the variables in the
5687 : * tableau have value "M - x" rather than "M + x".
5688 : */
5689 0 : static struct isl_sol *sol_pma_init(__isl_keep isl_basic_map *bmap,
5690 : __isl_take isl_basic_set *dom, int track_empty, int max)
5691 : {
5692 0 : struct isl_sol_pma *sol_pma = NULL;
5693 : isl_space *space;
5694 :
5695 0 : if (!bmap)
5696 0 : goto error;
5697 :
5698 0 : sol_pma = isl_calloc_type(bmap->ctx, struct isl_sol_pma);
5699 0 : if (!sol_pma)
5700 0 : goto error;
5701 :
5702 0 : sol_pma->sol.free = &sol_pma_free;
5703 0 : if (sol_init(&sol_pma->sol, bmap, dom, max) < 0)
5704 0 : goto error;
5705 0 : sol_pma->sol.add = &sol_pma_add_wrap;
5706 0 : sol_pma->sol.add_empty = track_empty ? &sol_pma_add_empty_wrap : NULL;
5707 0 : space = isl_space_copy(sol_pma->sol.space);
5708 0 : sol_pma->pma = isl_pw_multi_aff_empty(space);
5709 0 : if (!sol_pma->pma)
5710 0 : goto error;
5711 :
5712 0 : if (track_empty) {
5713 0 : sol_pma->empty = isl_set_alloc_space(isl_basic_set_get_space(dom),
5714 : 1, ISL_SET_DISJOINT);
5715 0 : if (!sol_pma->empty)
5716 0 : goto error;
5717 : }
5718 :
5719 0 : isl_basic_set_free(dom);
5720 0 : return &sol_pma->sol;
5721 : error:
5722 0 : isl_basic_set_free(dom);
5723 0 : sol_free(&sol_pma->sol);
5724 0 : return NULL;
5725 : }
5726 :
5727 : /* Base case of isl_tab_basic_map_partial_lexopt, after removing
5728 : * some obvious symmetries.
5729 : *
5730 : * We call basic_map_partial_lexopt_base_sol and extract the results.
5731 : */
5732 0 : static __isl_give isl_pw_multi_aff *basic_map_partial_lexopt_base_pw_multi_aff(
5733 : __isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom,
5734 : __isl_give isl_set **empty, int max)
5735 : {
5736 0 : isl_pw_multi_aff *result = NULL;
5737 : struct isl_sol *sol;
5738 : struct isl_sol_pma *sol_pma;
5739 :
5740 0 : sol = basic_map_partial_lexopt_base_sol(bmap, dom, empty, max,
5741 : &sol_pma_init);
5742 0 : if (!sol)
5743 0 : return NULL;
5744 0 : sol_pma = (struct isl_sol_pma *) sol;
5745 :
5746 0 : result = isl_pw_multi_aff_copy(sol_pma->pma);
5747 0 : if (empty)
5748 0 : *empty = isl_set_copy(sol_pma->empty);
5749 0 : sol_free(&sol_pma->sol);
5750 0 : return result;
5751 : }
5752 :
5753 : /* Given that the last input variable of "maff" represents the minimum
5754 : * of some bounds, check whether we need to plug in the expression
5755 : * of the minimum.
5756 : *
5757 : * In particular, check if the last input variable appears in any
5758 : * of the expressions in "maff".
5759 : */
5760 0 : static isl_bool need_substitution(__isl_keep isl_multi_aff *maff)
5761 : {
5762 : int i;
5763 : unsigned pos;
5764 :
5765 0 : pos = isl_multi_aff_dim(maff, isl_dim_in) - 1;
5766 :
5767 0 : for (i = 0; i < maff->n; ++i) {
5768 : isl_bool involves;
5769 :
5770 0 : involves = isl_aff_involves_dims(maff->u.p[i],
5771 : isl_dim_in, pos, 1);
5772 0 : if (involves < 0 || involves)
5773 0 : return involves;
5774 : }
5775 :
5776 0 : return isl_bool_false;
5777 : }
5778 :
5779 : /* Given a set of upper bounds on the last "input" variable m,
5780 : * construct a piecewise affine expression that selects
5781 : * the minimal upper bound to m, i.e.,
5782 : * divide the space into cells where one
5783 : * of the upper bounds is smaller than all the others and select
5784 : * this upper bound on that cell.
5785 : *
5786 : * In particular, if there are n bounds b_i, then the result
5787 : * consists of n cell, each one of the form
5788 : *
5789 : * b_i <= b_j for j > i
5790 : * b_i < b_j for j < i
5791 : *
5792 : * The affine expression on this cell is
5793 : *
5794 : * b_i
5795 : */
5796 0 : static __isl_give isl_pw_aff *set_minimum_pa(__isl_take isl_space *space,
5797 : __isl_take isl_mat *var)
5798 : {
5799 : int i;
5800 0 : isl_aff *aff = NULL;
5801 0 : isl_basic_set *bset = NULL;
5802 0 : isl_pw_aff *paff = NULL;
5803 : isl_space *pw_space;
5804 0 : isl_local_space *ls = NULL;
5805 :
5806 0 : if (!space || !var)
5807 : goto error;
5808 :
5809 0 : ls = isl_local_space_from_space(isl_space_copy(space));
5810 0 : pw_space = isl_space_copy(space);
5811 0 : pw_space = isl_space_from_domain(pw_space);
5812 0 : pw_space = isl_space_add_dims(pw_space, isl_dim_out, 1);
5813 0 : paff = isl_pw_aff_alloc_size(pw_space, var->n_row);
5814 :
5815 0 : for (i = 0; i < var->n_row; ++i) {
5816 : isl_pw_aff *paff_i;
5817 :
5818 0 : aff = isl_aff_alloc(isl_local_space_copy(ls));
5819 0 : bset = isl_basic_set_alloc_space(isl_space_copy(space), 0,
5820 0 : 0, var->n_row - 1);
5821 0 : if (!aff || !bset)
5822 : goto error;
5823 0 : isl_int_set_si(aff->v->el[0], 1);
5824 0 : isl_seq_cpy(aff->v->el + 1, var->row[i], var->n_col);
5825 0 : isl_int_set_si(aff->v->el[1 + var->n_col], 0);
5826 0 : bset = select_minimum(bset, var, i);
5827 0 : paff_i = isl_pw_aff_alloc(isl_set_from_basic_set(bset), aff);
5828 0 : paff = isl_pw_aff_add_disjoint(paff, paff_i);
5829 : }
5830 :
5831 0 : isl_local_space_free(ls);
5832 0 : isl_space_free(space);
5833 0 : isl_mat_free(var);
5834 0 : return paff;
5835 : error:
5836 0 : isl_aff_free(aff);
5837 0 : isl_basic_set_free(bset);
5838 0 : isl_pw_aff_free(paff);
5839 0 : isl_local_space_free(ls);
5840 0 : isl_space_free(space);
5841 0 : isl_mat_free(var);
5842 0 : return NULL;
5843 : }
5844 :
5845 : /* Given a piecewise multi-affine expression of which the last input variable
5846 : * is the minimum of the bounds in "cst", plug in the value of the minimum.
5847 : * This minimum expression is given in "min_expr_pa".
5848 : * The set "min_expr" contains the same information, but in the form of a set.
5849 : * The variable is subsequently projected out.
5850 : *
5851 : * The implementation is similar to those of "split" and "split_domain".
5852 : * If the variable appears in a given expression, then minimum expression
5853 : * is plugged in. Otherwise, if the variable appears in the constraints
5854 : * and a split is required, then the domain is split. Otherwise, no split
5855 : * is performed.
5856 : */
5857 0 : static __isl_give isl_pw_multi_aff *split_domain_pma(
5858 : __isl_take isl_pw_multi_aff *opt, __isl_take isl_pw_aff *min_expr_pa,
5859 : __isl_take isl_set *min_expr, __isl_take isl_mat *cst)
5860 : {
5861 : int n_in;
5862 : int i;
5863 : isl_space *space;
5864 : isl_pw_multi_aff *res;
5865 :
5866 0 : if (!opt || !min_expr || !cst)
5867 : goto error;
5868 :
5869 0 : n_in = isl_pw_multi_aff_dim(opt, isl_dim_in);
5870 0 : space = isl_pw_multi_aff_get_space(opt);
5871 0 : space = isl_space_drop_dims(space, isl_dim_in, n_in - 1, 1);
5872 0 : res = isl_pw_multi_aff_empty(space);
5873 :
5874 0 : for (i = 0; i < opt->n; ++i) {
5875 : isl_bool subs;
5876 : isl_pw_multi_aff *pma;
5877 :
5878 0 : pma = isl_pw_multi_aff_alloc(isl_set_copy(opt->p[i].set),
5879 : isl_multi_aff_copy(opt->p[i].maff));
5880 0 : subs = need_substitution(opt->p[i].maff);
5881 0 : if (subs < 0) {
5882 0 : pma = isl_pw_multi_aff_free(pma);
5883 0 : } else if (subs) {
5884 0 : pma = isl_pw_multi_aff_substitute(pma,
5885 0 : isl_dim_in, n_in - 1, min_expr_pa);
5886 : } else {
5887 : isl_bool split;
5888 0 : split = need_split_set(opt->p[i].set, cst);
5889 0 : if (split < 0)
5890 0 : pma = isl_pw_multi_aff_free(pma);
5891 0 : else if (split)
5892 0 : pma = isl_pw_multi_aff_intersect_domain(pma,
5893 : isl_set_copy(min_expr));
5894 : }
5895 0 : pma = isl_pw_multi_aff_project_out(pma,
5896 0 : isl_dim_in, n_in - 1, 1);
5897 :
5898 0 : res = isl_pw_multi_aff_add_disjoint(res, pma);
5899 : }
5900 :
5901 0 : isl_pw_multi_aff_free(opt);
5902 0 : isl_pw_aff_free(min_expr_pa);
5903 0 : isl_set_free(min_expr);
5904 0 : isl_mat_free(cst);
5905 0 : return res;
5906 : error:
5907 0 : isl_pw_multi_aff_free(opt);
5908 0 : isl_pw_aff_free(min_expr_pa);
5909 0 : isl_set_free(min_expr);
5910 0 : isl_mat_free(cst);
5911 0 : return NULL;
5912 : }
5913 :
5914 : static __isl_give isl_pw_multi_aff *basic_map_partial_lexopt_pw_multi_aff(
5915 : __isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom,
5916 : __isl_give isl_set **empty, int max);
5917 :
5918 : /* This function is called from basic_map_partial_lexopt_symm.
5919 : * The last variable of "bmap" and "dom" corresponds to the minimum
5920 : * of the bounds in "cst". "map_space" is the space of the original
5921 : * input relation (of basic_map_partial_lexopt_symm) and "set_space"
5922 : * is the space of the original domain.
5923 : *
5924 : * We recursively call basic_map_partial_lexopt and then plug in
5925 : * the definition of the minimum in the result.
5926 : */
5927 : static __isl_give isl_pw_multi_aff *
5928 0 : basic_map_partial_lexopt_symm_core_pw_multi_aff(
5929 : __isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom,
5930 : __isl_give isl_set **empty, int max, __isl_take isl_mat *cst,
5931 : __isl_take isl_space *map_space, __isl_take isl_space *set_space)
5932 : {
5933 : isl_pw_multi_aff *opt;
5934 : isl_pw_aff *min_expr_pa;
5935 : isl_set *min_expr;
5936 :
5937 0 : min_expr = set_minimum(isl_basic_set_get_space(dom), isl_mat_copy(cst));
5938 0 : min_expr_pa = set_minimum_pa(isl_basic_set_get_space(dom),
5939 : isl_mat_copy(cst));
5940 :
5941 0 : opt = basic_map_partial_lexopt_pw_multi_aff(bmap, dom, empty, max);
5942 :
5943 0 : if (empty) {
5944 0 : *empty = split(*empty,
5945 : isl_set_copy(min_expr), isl_mat_copy(cst));
5946 0 : *empty = isl_set_reset_space(*empty, set_space);
5947 : }
5948 :
5949 0 : opt = split_domain_pma(opt, min_expr_pa, min_expr, cst);
5950 0 : opt = isl_pw_multi_aff_reset_space(opt, map_space);
5951 :
5952 0 : return opt;
5953 : }
5954 :
5955 : #undef TYPE
5956 : #define TYPE isl_pw_multi_aff
5957 : #undef SUFFIX
5958 : #define SUFFIX _pw_multi_aff
5959 : #include "isl_tab_lexopt_templ.c"
|
