add isl_pw_qpolynomial_bound_range
[isl.git] / isl_tab_pip.c
blob36c3d02e2f3b81f9feb63d61889120f2e36d8e4c
1 /*
2 * Copyright 2008-2009 Katholieke Universiteit Leuven
4 * Use of this software is governed by the GNU LGPLv2.1 license
6 * Written by Sven Verdoolaege, K.U.Leuven, Departement
7 * Computerwetenschappen, Celestijnenlaan 200A, B-3001 Leuven, Belgium
8 */
10 #include "isl_map_private.h"
11 #include "isl_seq.h"
12 #include "isl_tab.h"
13 #include "isl_sample.h"
16 * The implementation of parametric integer linear programming in this file
17 * was inspired by the paper "Parametric Integer Programming" and the
18 * report "Solving systems of affine (in)equalities" by Paul Feautrier
19 * (and others).
21 * The strategy used for obtaining a feasible solution is different
22 * from the one used in isl_tab.c. In particular, in isl_tab.c,
23 * upon finding a constraint that is not yet satisfied, we pivot
24 * in a row that increases the constant term of row holding the
25 * constraint, making sure the sample solution remains feasible
26 * for all the constraints it already satisfied.
27 * Here, we always pivot in the row holding the constraint,
28 * choosing a column that induces the lexicographically smallest
29 * increment to the sample solution.
31 * By starting out from a sample value that is lexicographically
32 * smaller than any integer point in the problem space, the first
33 * feasible integer sample point we find will also be the lexicographically
34 * smallest. If all variables can be assumed to be non-negative,
35 * then the initial sample value may be chosen equal to zero.
36 * However, we will not make this assumption. Instead, we apply
37 * the "big parameter" trick. Any variable x is then not directly
38 * used in the tableau, but instead it its represented by another
39 * variable x' = M + x, where M is an arbitrarily large (positive)
40 * value. x' is therefore always non-negative, whatever the value of x.
41 * Taking as initial smaple value x' = 0 corresponds to x = -M,
42 * which is always smaller than any possible value of x.
44 * The big parameter trick is used in the main tableau and
45 * also in the context tableau if isl_context_lex is used.
46 * In this case, each tableaus has its own big parameter.
47 * Before doing any real work, we check if all the parameters
48 * happen to be non-negative. If so, we drop the column corresponding
49 * to M from the initial context tableau.
50 * If isl_context_gbr is used, then the big parameter trick is only
51 * used in the main tableau.
54 struct isl_context;
55 struct isl_context_op {
56 /* detect nonnegative parameters in context and mark them in tab */
57 struct isl_tab *(*detect_nonnegative_parameters)(
58 struct isl_context *context, struct isl_tab *tab);
59 /* return temporary reference to basic set representation of context */
60 struct isl_basic_set *(*peek_basic_set)(struct isl_context *context);
61 /* return temporary reference to tableau representation of context */
62 struct isl_tab *(*peek_tab)(struct isl_context *context);
63 /* add equality; check is 1 if eq may not be valid;
64 * update is 1 if we may want to call ineq_sign on context later.
66 void (*add_eq)(struct isl_context *context, isl_int *eq,
67 int check, int update);
68 /* add inequality; check is 1 if ineq may not be valid;
69 * update is 1 if we may want to call ineq_sign on context later.
71 void (*add_ineq)(struct isl_context *context, isl_int *ineq,
72 int check, int update);
73 /* check sign of ineq based on previous information.
74 * strict is 1 if saturation should be treated as a positive sign.
76 enum isl_tab_row_sign (*ineq_sign)(struct isl_context *context,
77 isl_int *ineq, int strict);
78 /* check if inequality maintains feasibility */
79 int (*test_ineq)(struct isl_context *context, isl_int *ineq);
80 /* return index of a div that corresponds to "div" */
81 int (*get_div)(struct isl_context *context, struct isl_tab *tab,
82 struct isl_vec *div);
83 /* add div "div" to context and return non-negativity */
84 int (*add_div)(struct isl_context *context, struct isl_vec *div);
85 int (*detect_equalities)(struct isl_context *context,
86 struct isl_tab *tab);
87 /* return row index of "best" split */
88 int (*best_split)(struct isl_context *context, struct isl_tab *tab);
89 /* check if context has already been determined to be empty */
90 int (*is_empty)(struct isl_context *context);
91 /* check if context is still usable */
92 int (*is_ok)(struct isl_context *context);
93 /* save a copy/snapshot of context */
94 void *(*save)(struct isl_context *context);
95 /* restore saved context */
96 void (*restore)(struct isl_context *context, void *);
97 /* invalidate context */
98 void (*invalidate)(struct isl_context *context);
99 /* free context */
100 void (*free)(struct isl_context *context);
103 struct isl_context {
104 struct isl_context_op *op;
107 struct isl_context_lex {
108 struct isl_context context;
109 struct isl_tab *tab;
112 struct isl_partial_sol {
113 int level;
114 struct isl_basic_set *dom;
115 struct isl_mat *M;
117 struct isl_partial_sol *next;
120 struct isl_sol;
121 struct isl_sol_callback {
122 struct isl_tab_callback callback;
123 struct isl_sol *sol;
126 /* isl_sol is an interface for constructing a solution to
127 * a parametric integer linear programming problem.
128 * Every time the algorithm reaches a state where a solution
129 * can be read off from the tableau (including cases where the tableau
130 * is empty), the function "add" is called on the isl_sol passed
131 * to find_solutions_main.
133 * The context tableau is owned by isl_sol and is updated incrementally.
135 * There are currently two implementations of this interface,
136 * isl_sol_map, which simply collects the solutions in an isl_map
137 * and (optionally) the parts of the context where there is no solution
138 * in an isl_set, and
139 * isl_sol_for, which calls a user-defined function for each part of
140 * the solution.
142 struct isl_sol {
143 int error;
144 int rational;
145 int level;
146 int max;
147 int n_out;
148 struct isl_context *context;
149 struct isl_partial_sol *partial;
150 void (*add)(struct isl_sol *sol,
151 struct isl_basic_set *dom, struct isl_mat *M);
152 void (*add_empty)(struct isl_sol *sol, struct isl_basic_set *bset);
153 void (*free)(struct isl_sol *sol);
154 struct isl_sol_callback dec_level;
157 static void sol_free(struct isl_sol *sol)
159 struct isl_partial_sol *partial, *next;
160 if (!sol)
161 return;
162 for (partial = sol->partial; partial; partial = next) {
163 next = partial->next;
164 isl_basic_set_free(partial->dom);
165 isl_mat_free(partial->M);
166 free(partial);
168 sol->free(sol);
171 /* Push a partial solution represented by a domain and mapping M
172 * onto the stack of partial solutions.
174 static void sol_push_sol(struct isl_sol *sol,
175 struct isl_basic_set *dom, struct isl_mat *M)
177 struct isl_partial_sol *partial;
179 if (sol->error || !dom)
180 goto error;
182 partial = isl_alloc_type(dom->ctx, struct isl_partial_sol);
183 if (!partial)
184 goto error;
186 partial->level = sol->level;
187 partial->dom = dom;
188 partial->M = M;
189 partial->next = sol->partial;
191 sol->partial = partial;
193 return;
194 error:
195 isl_basic_set_free(dom);
196 sol->error = 1;
199 /* Pop one partial solution from the partial solution stack and
200 * pass it on to sol->add or sol->add_empty.
202 static void sol_pop_one(struct isl_sol *sol)
204 struct isl_partial_sol *partial;
206 partial = sol->partial;
207 sol->partial = partial->next;
209 if (partial->M)
210 sol->add(sol, partial->dom, partial->M);
211 else
212 sol->add_empty(sol, partial->dom);
213 free(partial);
216 /* Return a fresh copy of the domain represented by the context tableau.
218 static struct isl_basic_set *sol_domain(struct isl_sol *sol)
220 struct isl_basic_set *bset;
222 if (sol->error)
223 return NULL;
225 bset = isl_basic_set_dup(sol->context->op->peek_basic_set(sol->context));
226 bset = isl_basic_set_update_from_tab(bset,
227 sol->context->op->peek_tab(sol->context));
229 return bset;
232 /* Check whether two partial solutions have the same mapping, where n_div
233 * is the number of divs that the two partial solutions have in common.
235 static int same_solution(struct isl_partial_sol *s1, struct isl_partial_sol *s2,
236 unsigned n_div)
238 int i;
239 unsigned dim;
241 if (!s1->M != !s2->M)
242 return 0;
243 if (!s1->M)
244 return 1;
246 dim = isl_basic_set_total_dim(s1->dom) - s1->dom->n_div;
248 for (i = 0; i < s1->M->n_row; ++i) {
249 if (isl_seq_first_non_zero(s1->M->row[i]+1+dim+n_div,
250 s1->M->n_col-1-dim-n_div) != -1)
251 return 0;
252 if (isl_seq_first_non_zero(s2->M->row[i]+1+dim+n_div,
253 s2->M->n_col-1-dim-n_div) != -1)
254 return 0;
255 if (!isl_seq_eq(s1->M->row[i], s2->M->row[i], 1+dim+n_div))
256 return 0;
258 return 1;
261 /* Pop all solutions from the partial solution stack that were pushed onto
262 * the stack at levels that are deeper than the current level.
263 * If the two topmost elements on the stack have the same level
264 * and represent the same solution, then their domains are combined.
265 * This combined domain is the same as the current context domain
266 * as sol_pop is called each time we move back to a higher level.
268 static void sol_pop(struct isl_sol *sol)
270 struct isl_partial_sol *partial;
271 unsigned n_div;
273 if (sol->error)
274 return;
276 if (sol->level == 0) {
277 for (partial = sol->partial; partial; partial = sol->partial)
278 sol_pop_one(sol);
279 return;
282 partial = sol->partial;
283 if (!partial)
284 return;
286 if (partial->level <= sol->level)
287 return;
289 if (partial->next && partial->next->level == partial->level) {
290 n_div = isl_basic_set_dim(
291 sol->context->op->peek_basic_set(sol->context),
292 isl_dim_div);
294 if (!same_solution(partial, partial->next, n_div)) {
295 sol_pop_one(sol);
296 sol_pop_one(sol);
297 } else {
298 struct isl_basic_set *bset;
300 bset = sol_domain(sol);
302 isl_basic_set_free(partial->next->dom);
303 partial->next->dom = bset;
304 partial->next->level = sol->level;
306 sol->partial = partial->next;
307 isl_basic_set_free(partial->dom);
308 isl_mat_free(partial->M);
309 free(partial);
311 } else
312 sol_pop_one(sol);
315 static void sol_dec_level(struct isl_sol *sol)
317 if (sol->error)
318 return;
320 sol->level--;
322 sol_pop(sol);
325 static int sol_dec_level_wrap(struct isl_tab_callback *cb)
327 struct isl_sol_callback *callback = (struct isl_sol_callback *)cb;
329 sol_dec_level(callback->sol);
331 return callback->sol->error ? -1 : 0;
334 /* Move down to next level and push callback onto context tableau
335 * to decrease the level again when it gets rolled back across
336 * the current state. That is, dec_level will be called with
337 * the context tableau in the same state as it is when inc_level
338 * is called.
340 static void sol_inc_level(struct isl_sol *sol)
342 struct isl_tab *tab;
344 if (sol->error)
345 return;
347 sol->level++;
348 tab = sol->context->op->peek_tab(sol->context);
349 if (isl_tab_push_callback(tab, &sol->dec_level.callback) < 0)
350 sol->error = 1;
353 static void scale_rows(struct isl_mat *mat, isl_int m, int n_row)
355 int i;
357 if (isl_int_is_one(m))
358 return;
360 for (i = 0; i < n_row; ++i)
361 isl_seq_scale(mat->row[i], mat->row[i], m, mat->n_col);
364 /* Add the solution identified by the tableau and the context tableau.
366 * The layout of the variables is as follows.
367 * tab->n_var is equal to the total number of variables in the input
368 * map (including divs that were copied from the context)
369 * + the number of extra divs constructed
370 * Of these, the first tab->n_param and the last tab->n_div variables
371 * correspond to the variables in the context, i.e.,
372 * tab->n_param + tab->n_div = context_tab->n_var
373 * tab->n_param is equal to the number of parameters and input
374 * dimensions in the input map
375 * tab->n_div is equal to the number of divs in the context
377 * If there is no solution, then call add_empty with a basic set
378 * that corresponds to the context tableau. (If add_empty is NULL,
379 * then do nothing).
381 * If there is a solution, then first construct a matrix that maps
382 * all dimensions of the context to the output variables, i.e.,
383 * the output dimensions in the input map.
384 * The divs in the input map (if any) that do not correspond to any
385 * div in the context do not appear in the solution.
386 * The algorithm will make sure that they have an integer value,
387 * but these values themselves are of no interest.
388 * We have to be careful not to drop or rearrange any divs in the
389 * context because that would change the meaning of the matrix.
391 * To extract the value of the output variables, it should be noted
392 * that we always use a big parameter M in the main tableau and so
393 * the variable stored in this tableau is not an output variable x itself, but
394 * x' = M + x (in case of minimization)
395 * or
396 * x' = M - x (in case of maximization)
397 * If x' appears in a column, then its optimal value is zero,
398 * which means that the optimal value of x is an unbounded number
399 * (-M for minimization and M for maximization).
400 * We currently assume that the output dimensions in the original map
401 * are bounded, so this cannot occur.
402 * Similarly, when x' appears in a row, then the coefficient of M in that
403 * row is necessarily 1.
404 * If the row in the tableau represents
405 * d x' = c + d M + e(y)
406 * then, in case of minimization, the corresponding row in the matrix
407 * will be
408 * a c + a e(y)
409 * with a d = m, the (updated) common denominator of the matrix.
410 * In case of maximization, the row will be
411 * -a c - a e(y)
413 static void sol_add(struct isl_sol *sol, struct isl_tab *tab)
415 struct isl_basic_set *bset = NULL;
416 struct isl_mat *mat = NULL;
417 unsigned off;
418 int row, i;
419 isl_int m;
421 if (sol->error || !tab)
422 goto error;
424 if (tab->empty && !sol->add_empty)
425 return;
427 bset = sol_domain(sol);
429 if (tab->empty) {
430 sol_push_sol(sol, bset, NULL);
431 return;
434 off = 2 + tab->M;
436 mat = isl_mat_alloc(tab->mat->ctx, 1 + sol->n_out,
437 1 + tab->n_param + tab->n_div);
438 if (!mat)
439 goto error;
441 isl_int_init(m);
443 isl_seq_clr(mat->row[0] + 1, mat->n_col - 1);
444 isl_int_set_si(mat->row[0][0], 1);
445 for (row = 0; row < sol->n_out; ++row) {
446 int i = tab->n_param + row;
447 int r, j;
449 isl_seq_clr(mat->row[1 + row], mat->n_col);
450 if (!tab->var[i].is_row) {
451 /* no unbounded */
452 isl_assert(mat->ctx, !tab->M, goto error2);
453 continue;
456 r = tab->var[i].index;
457 /* no unbounded */
458 if (tab->M)
459 isl_assert(mat->ctx, isl_int_eq(tab->mat->row[r][2],
460 tab->mat->row[r][0]),
461 goto error2);
462 isl_int_gcd(m, mat->row[0][0], tab->mat->row[r][0]);
463 isl_int_divexact(m, tab->mat->row[r][0], m);
464 scale_rows(mat, m, 1 + row);
465 isl_int_divexact(m, mat->row[0][0], tab->mat->row[r][0]);
466 isl_int_mul(mat->row[1 + row][0], m, tab->mat->row[r][1]);
467 for (j = 0; j < tab->n_param; ++j) {
468 int col;
469 if (tab->var[j].is_row)
470 continue;
471 col = tab->var[j].index;
472 isl_int_mul(mat->row[1 + row][1 + j], m,
473 tab->mat->row[r][off + col]);
475 for (j = 0; j < tab->n_div; ++j) {
476 int col;
477 if (tab->var[tab->n_var - tab->n_div+j].is_row)
478 continue;
479 col = tab->var[tab->n_var - tab->n_div+j].index;
480 isl_int_mul(mat->row[1 + row][1 + tab->n_param + j], m,
481 tab->mat->row[r][off + col]);
483 if (sol->max)
484 isl_seq_neg(mat->row[1 + row], mat->row[1 + row],
485 mat->n_col);
488 isl_int_clear(m);
490 sol_push_sol(sol, bset, mat);
491 return;
492 error2:
493 isl_int_clear(m);
494 error:
495 isl_basic_set_free(bset);
496 isl_mat_free(mat);
497 sol_free(sol);
500 struct isl_sol_map {
501 struct isl_sol sol;
502 struct isl_map *map;
503 struct isl_set *empty;
506 static void sol_map_free(struct isl_sol_map *sol_map)
508 if (sol_map->sol.context)
509 sol_map->sol.context->op->free(sol_map->sol.context);
510 isl_map_free(sol_map->map);
511 isl_set_free(sol_map->empty);
512 free(sol_map);
515 static void sol_map_free_wrap(struct isl_sol *sol)
517 sol_map_free((struct isl_sol_map *)sol);
520 /* This function is called for parts of the context where there is
521 * no solution, with "bset" corresponding to the context tableau.
522 * Simply add the basic set to the set "empty".
524 static void sol_map_add_empty(struct isl_sol_map *sol,
525 struct isl_basic_set *bset)
527 if (!bset)
528 goto error;
529 isl_assert(bset->ctx, sol->empty, goto error);
531 sol->empty = isl_set_grow(sol->empty, 1);
532 bset = isl_basic_set_simplify(bset);
533 bset = isl_basic_set_finalize(bset);
534 sol->empty = isl_set_add_basic_set(sol->empty, isl_basic_set_copy(bset));
535 if (!sol->empty)
536 goto error;
537 isl_basic_set_free(bset);
538 return;
539 error:
540 isl_basic_set_free(bset);
541 sol->sol.error = 1;
544 static void sol_map_add_empty_wrap(struct isl_sol *sol,
545 struct isl_basic_set *bset)
547 sol_map_add_empty((struct isl_sol_map *)sol, bset);
550 /* Add bset to sol's empty, but only if we are actually collecting
551 * the empty set.
553 static void sol_map_add_empty_if_needed(struct isl_sol_map *sol,
554 struct isl_basic_set *bset)
556 if (sol->empty)
557 sol_map_add_empty(sol, bset);
558 else
559 isl_basic_set_free(bset);
562 /* Given a basic map "dom" that represents the context and an affine
563 * matrix "M" that maps the dimensions of the context to the
564 * output variables, construct a basic map with the same parameters
565 * and divs as the context, the dimensions of the context as input
566 * dimensions and a number of output dimensions that is equal to
567 * the number of output dimensions in the input map.
569 * The constraints and divs of the context are simply copied
570 * from "dom". For each row
571 * x = c + e(y)
572 * an equality
573 * c + e(y) - d x = 0
574 * is added, with d the common denominator of M.
576 static void sol_map_add(struct isl_sol_map *sol,
577 struct isl_basic_set *dom, struct isl_mat *M)
579 int i;
580 struct isl_basic_map *bmap = NULL;
581 isl_basic_set *context_bset;
582 unsigned n_eq;
583 unsigned n_ineq;
584 unsigned nparam;
585 unsigned total;
586 unsigned n_div;
587 unsigned n_out;
589 if (sol->sol.error || !dom || !M)
590 goto error;
592 n_out = sol->sol.n_out;
593 n_eq = dom->n_eq + n_out;
594 n_ineq = dom->n_ineq;
595 n_div = dom->n_div;
596 nparam = isl_basic_set_total_dim(dom) - n_div;
597 total = isl_map_dim(sol->map, isl_dim_all);
598 bmap = isl_basic_map_alloc_dim(isl_map_get_dim(sol->map),
599 n_div, n_eq, 2 * n_div + n_ineq);
600 if (!bmap)
601 goto error;
602 if (sol->sol.rational)
603 ISL_F_SET(bmap, ISL_BASIC_MAP_RATIONAL);
604 for (i = 0; i < dom->n_div; ++i) {
605 int k = isl_basic_map_alloc_div(bmap);
606 if (k < 0)
607 goto error;
608 isl_seq_cpy(bmap->div[k], dom->div[i], 1 + 1 + nparam);
609 isl_seq_clr(bmap->div[k] + 1 + 1 + nparam, total - nparam);
610 isl_seq_cpy(bmap->div[k] + 1 + 1 + total,
611 dom->div[i] + 1 + 1 + nparam, i);
613 for (i = 0; i < dom->n_eq; ++i) {
614 int k = isl_basic_map_alloc_equality(bmap);
615 if (k < 0)
616 goto error;
617 isl_seq_cpy(bmap->eq[k], dom->eq[i], 1 + nparam);
618 isl_seq_clr(bmap->eq[k] + 1 + nparam, total - nparam);
619 isl_seq_cpy(bmap->eq[k] + 1 + total,
620 dom->eq[i] + 1 + nparam, n_div);
622 for (i = 0; i < dom->n_ineq; ++i) {
623 int k = isl_basic_map_alloc_inequality(bmap);
624 if (k < 0)
625 goto error;
626 isl_seq_cpy(bmap->ineq[k], dom->ineq[i], 1 + nparam);
627 isl_seq_clr(bmap->ineq[k] + 1 + nparam, total - nparam);
628 isl_seq_cpy(bmap->ineq[k] + 1 + total,
629 dom->ineq[i] + 1 + nparam, n_div);
631 for (i = 0; i < M->n_row - 1; ++i) {
632 int k = isl_basic_map_alloc_equality(bmap);
633 if (k < 0)
634 goto error;
635 isl_seq_cpy(bmap->eq[k], M->row[1 + i], 1 + nparam);
636 isl_seq_clr(bmap->eq[k] + 1 + nparam, n_out);
637 isl_int_neg(bmap->eq[k][1 + nparam + i], M->row[0][0]);
638 isl_seq_cpy(bmap->eq[k] + 1 + nparam + n_out,
639 M->row[1 + i] + 1 + nparam, n_div);
641 bmap = isl_basic_map_simplify(bmap);
642 bmap = isl_basic_map_finalize(bmap);
643 sol->map = isl_map_grow(sol->map, 1);
644 sol->map = isl_map_add_basic_map(sol->map, bmap);
645 if (!sol->map)
646 goto error;
647 isl_basic_set_free(dom);
648 isl_mat_free(M);
649 return;
650 error:
651 isl_basic_set_free(dom);
652 isl_mat_free(M);
653 isl_basic_map_free(bmap);
654 sol->sol.error = 1;
657 static void sol_map_add_wrap(struct isl_sol *sol,
658 struct isl_basic_set *dom, struct isl_mat *M)
660 sol_map_add((struct isl_sol_map *)sol, dom, M);
664 /* Store the "parametric constant" of row "row" of tableau "tab" in "line",
665 * i.e., the constant term and the coefficients of all variables that
666 * appear in the context tableau.
667 * Note that the coefficient of the big parameter M is NOT copied.
668 * The context tableau may not have a big parameter and even when it
669 * does, it is a different big parameter.
671 static void get_row_parameter_line(struct isl_tab *tab, int row, isl_int *line)
673 int i;
674 unsigned off = 2 + tab->M;
676 isl_int_set(line[0], tab->mat->row[row][1]);
677 for (i = 0; i < tab->n_param; ++i) {
678 if (tab->var[i].is_row)
679 isl_int_set_si(line[1 + i], 0);
680 else {
681 int col = tab->var[i].index;
682 isl_int_set(line[1 + i], tab->mat->row[row][off + col]);
685 for (i = 0; i < tab->n_div; ++i) {
686 if (tab->var[tab->n_var - tab->n_div + i].is_row)
687 isl_int_set_si(line[1 + tab->n_param + i], 0);
688 else {
689 int col = tab->var[tab->n_var - tab->n_div + i].index;
690 isl_int_set(line[1 + tab->n_param + i],
691 tab->mat->row[row][off + col]);
696 /* Check if rows "row1" and "row2" have identical "parametric constants",
697 * as explained above.
698 * In this case, we also insist that the coefficients of the big parameter
699 * be the same as the values of the constants will only be the same
700 * if these coefficients are also the same.
702 static int identical_parameter_line(struct isl_tab *tab, int row1, int row2)
704 int i;
705 unsigned off = 2 + tab->M;
707 if (isl_int_ne(tab->mat->row[row1][1], tab->mat->row[row2][1]))
708 return 0;
710 if (tab->M && isl_int_ne(tab->mat->row[row1][2],
711 tab->mat->row[row2][2]))
712 return 0;
714 for (i = 0; i < tab->n_param + tab->n_div; ++i) {
715 int pos = i < tab->n_param ? i :
716 tab->n_var - tab->n_div + i - tab->n_param;
717 int col;
719 if (tab->var[pos].is_row)
720 continue;
721 col = tab->var[pos].index;
722 if (isl_int_ne(tab->mat->row[row1][off + col],
723 tab->mat->row[row2][off + col]))
724 return 0;
726 return 1;
729 /* Return an inequality that expresses that the "parametric constant"
730 * should be non-negative.
731 * This function is only called when the coefficient of the big parameter
732 * is equal to zero.
734 static struct isl_vec *get_row_parameter_ineq(struct isl_tab *tab, int row)
736 struct isl_vec *ineq;
738 ineq = isl_vec_alloc(tab->mat->ctx, 1 + tab->n_param + tab->n_div);
739 if (!ineq)
740 return NULL;
742 get_row_parameter_line(tab, row, ineq->el);
743 if (ineq)
744 ineq = isl_vec_normalize(ineq);
746 return ineq;
749 /* Return a integer division for use in a parametric cut based on the given row.
750 * In particular, let the parametric constant of the row be
752 * \sum_i a_i y_i
754 * where y_0 = 1, but none of the y_i corresponds to the big parameter M.
755 * The div returned is equal to
757 * floor(\sum_i {-a_i} y_i) = floor((\sum_i (-a_i mod d) y_i)/d)
759 static struct isl_vec *get_row_parameter_div(struct isl_tab *tab, int row)
761 struct isl_vec *div;
763 div = isl_vec_alloc(tab->mat->ctx, 1 + 1 + tab->n_param + tab->n_div);
764 if (!div)
765 return NULL;
767 isl_int_set(div->el[0], tab->mat->row[row][0]);
768 get_row_parameter_line(tab, row, div->el + 1);
769 div = isl_vec_normalize(div);
770 isl_seq_neg(div->el + 1, div->el + 1, div->size - 1);
771 isl_seq_fdiv_r(div->el + 1, div->el + 1, div->el[0], div->size - 1);
773 return div;
776 /* Return a integer division for use in transferring an integrality constraint
777 * to the context.
778 * In particular, let the parametric constant of the row be
780 * \sum_i a_i y_i
782 * where y_0 = 1, but none of the y_i corresponds to the big parameter M.
783 * The the returned div is equal to
785 * floor(\sum_i {a_i} y_i) = floor((\sum_i (a_i mod d) y_i)/d)
787 static struct isl_vec *get_row_split_div(struct isl_tab *tab, int row)
789 struct isl_vec *div;
791 div = isl_vec_alloc(tab->mat->ctx, 1 + 1 + tab->n_param + tab->n_div);
792 if (!div)
793 return NULL;
795 isl_int_set(div->el[0], tab->mat->row[row][0]);
796 get_row_parameter_line(tab, row, div->el + 1);
797 div = isl_vec_normalize(div);
798 isl_seq_fdiv_r(div->el + 1, div->el + 1, div->el[0], div->size - 1);
800 return div;
803 /* Construct and return an inequality that expresses an upper bound
804 * on the given div.
805 * In particular, if the div is given by
807 * d = floor(e/m)
809 * then the inequality expresses
811 * m d <= e
813 static struct isl_vec *ineq_for_div(struct isl_basic_set *bset, unsigned div)
815 unsigned total;
816 unsigned div_pos;
817 struct isl_vec *ineq;
819 if (!bset)
820 return NULL;
822 total = isl_basic_set_total_dim(bset);
823 div_pos = 1 + total - bset->n_div + div;
825 ineq = isl_vec_alloc(bset->ctx, 1 + total);
826 if (!ineq)
827 return NULL;
829 isl_seq_cpy(ineq->el, bset->div[div] + 1, 1 + total);
830 isl_int_neg(ineq->el[div_pos], bset->div[div][0]);
831 return ineq;
834 /* Given a row in the tableau and a div that was created
835 * using get_row_split_div and that been constrained to equality, i.e.,
837 * d = floor(\sum_i {a_i} y_i) = \sum_i {a_i} y_i
839 * replace the expression "\sum_i {a_i} y_i" in the row by d,
840 * i.e., we subtract "\sum_i {a_i} y_i" and add 1 d.
841 * The coefficients of the non-parameters in the tableau have been
842 * verified to be integral. We can therefore simply replace coefficient b
843 * by floor(b). For the coefficients of the parameters we have
844 * floor(a_i) = a_i - {a_i}, while for the other coefficients, we have
845 * floor(b) = b.
847 static struct isl_tab *set_row_cst_to_div(struct isl_tab *tab, int row, int div)
849 isl_seq_fdiv_q(tab->mat->row[row] + 1, tab->mat->row[row] + 1,
850 tab->mat->row[row][0], 1 + tab->M + tab->n_col);
852 isl_int_set_si(tab->mat->row[row][0], 1);
854 if (tab->var[tab->n_var - tab->n_div + div].is_row) {
855 int drow = tab->var[tab->n_var - tab->n_div + div].index;
857 isl_assert(tab->mat->ctx,
858 isl_int_is_one(tab->mat->row[drow][0]), goto error);
859 isl_seq_combine(tab->mat->row[row] + 1,
860 tab->mat->ctx->one, tab->mat->row[row] + 1,
861 tab->mat->ctx->one, tab->mat->row[drow] + 1,
862 1 + tab->M + tab->n_col);
863 } else {
864 int dcol = tab->var[tab->n_var - tab->n_div + div].index;
866 isl_int_set_si(tab->mat->row[row][2 + tab->M + dcol], 1);
869 return tab;
870 error:
871 isl_tab_free(tab);
872 return NULL;
875 /* Check if the (parametric) constant of the given row is obviously
876 * negative, meaning that we don't need to consult the context tableau.
877 * If there is a big parameter and its coefficient is non-zero,
878 * then this coefficient determines the outcome.
879 * Otherwise, we check whether the constant is negative and
880 * all non-zero coefficients of parameters are negative and
881 * belong to non-negative parameters.
883 static int is_obviously_neg(struct isl_tab *tab, int row)
885 int i;
886 int col;
887 unsigned off = 2 + tab->M;
889 if (tab->M) {
890 if (isl_int_is_pos(tab->mat->row[row][2]))
891 return 0;
892 if (isl_int_is_neg(tab->mat->row[row][2]))
893 return 1;
896 if (isl_int_is_nonneg(tab->mat->row[row][1]))
897 return 0;
898 for (i = 0; i < tab->n_param; ++i) {
899 /* Eliminated parameter */
900 if (tab->var[i].is_row)
901 continue;
902 col = tab->var[i].index;
903 if (isl_int_is_zero(tab->mat->row[row][off + col]))
904 continue;
905 if (!tab->var[i].is_nonneg)
906 return 0;
907 if (isl_int_is_pos(tab->mat->row[row][off + col]))
908 return 0;
910 for (i = 0; i < tab->n_div; ++i) {
911 if (tab->var[tab->n_var - tab->n_div + i].is_row)
912 continue;
913 col = tab->var[tab->n_var - tab->n_div + i].index;
914 if (isl_int_is_zero(tab->mat->row[row][off + col]))
915 continue;
916 if (!tab->var[tab->n_var - tab->n_div + i].is_nonneg)
917 return 0;
918 if (isl_int_is_pos(tab->mat->row[row][off + col]))
919 return 0;
921 return 1;
924 /* Check if the (parametric) constant of the given row is obviously
925 * non-negative, meaning that we don't need to consult the context tableau.
926 * If there is a big parameter and its coefficient is non-zero,
927 * then this coefficient determines the outcome.
928 * Otherwise, we check whether the constant is non-negative and
929 * all non-zero coefficients of parameters are positive and
930 * belong to non-negative parameters.
932 static int is_obviously_nonneg(struct isl_tab *tab, int row)
934 int i;
935 int col;
936 unsigned off = 2 + tab->M;
938 if (tab->M) {
939 if (isl_int_is_pos(tab->mat->row[row][2]))
940 return 1;
941 if (isl_int_is_neg(tab->mat->row[row][2]))
942 return 0;
945 if (isl_int_is_neg(tab->mat->row[row][1]))
946 return 0;
947 for (i = 0; i < tab->n_param; ++i) {
948 /* Eliminated parameter */
949 if (tab->var[i].is_row)
950 continue;
951 col = tab->var[i].index;
952 if (isl_int_is_zero(tab->mat->row[row][off + col]))
953 continue;
954 if (!tab->var[i].is_nonneg)
955 return 0;
956 if (isl_int_is_neg(tab->mat->row[row][off + col]))
957 return 0;
959 for (i = 0; i < tab->n_div; ++i) {
960 if (tab->var[tab->n_var - tab->n_div + i].is_row)
961 continue;
962 col = tab->var[tab->n_var - tab->n_div + i].index;
963 if (isl_int_is_zero(tab->mat->row[row][off + col]))
964 continue;
965 if (!tab->var[tab->n_var - tab->n_div + i].is_nonneg)
966 return 0;
967 if (isl_int_is_neg(tab->mat->row[row][off + col]))
968 return 0;
970 return 1;
973 /* Given a row r and two columns, return the column that would
974 * lead to the lexicographically smallest increment in the sample
975 * solution when leaving the basis in favor of the row.
976 * Pivoting with column c will increment the sample value by a non-negative
977 * constant times a_{V,c}/a_{r,c}, with a_{V,c} the elements of column c
978 * corresponding to the non-parametric variables.
979 * If variable v appears in a column c_v, the a_{v,c} = 1 iff c = c_v,
980 * with all other entries in this virtual row equal to zero.
981 * If variable v appears in a row, then a_{v,c} is the element in column c
982 * of that row.
984 * Let v be the first variable with a_{v,c1}/a_{r,c1} != a_{v,c2}/a_{r,c2}.
985 * Then if a_{v,c1}/a_{r,c1} < a_{v,c2}/a_{r,c2}, i.e.,
986 * a_{v,c2} a_{r,c1} - a_{v,c1} a_{r,c2} > 0, c1 results in the minimal
987 * increment. Otherwise, it's c2.
989 static int lexmin_col_pair(struct isl_tab *tab,
990 int row, int col1, int col2, isl_int tmp)
992 int i;
993 isl_int *tr;
995 tr = tab->mat->row[row] + 2 + tab->M;
997 for (i = tab->n_param; i < tab->n_var - tab->n_div; ++i) {
998 int s1, s2;
999 isl_int *r;
1001 if (!tab->var[i].is_row) {
1002 if (tab->var[i].index == col1)
1003 return col2;
1004 if (tab->var[i].index == col2)
1005 return col1;
1006 continue;
1009 if (tab->var[i].index == row)
1010 continue;
1012 r = tab->mat->row[tab->var[i].index] + 2 + tab->M;
1013 s1 = isl_int_sgn(r[col1]);
1014 s2 = isl_int_sgn(r[col2]);
1015 if (s1 == 0 && s2 == 0)
1016 continue;
1017 if (s1 < s2)
1018 return col1;
1019 if (s2 < s1)
1020 return col2;
1022 isl_int_mul(tmp, r[col2], tr[col1]);
1023 isl_int_submul(tmp, r[col1], tr[col2]);
1024 if (isl_int_is_pos(tmp))
1025 return col1;
1026 if (isl_int_is_neg(tmp))
1027 return col2;
1029 return -1;
1032 /* Given a row in the tableau, find and return the column that would
1033 * result in the lexicographically smallest, but positive, increment
1034 * in the sample point.
1035 * If there is no such column, then return tab->n_col.
1036 * If anything goes wrong, return -1.
1038 static int lexmin_pivot_col(struct isl_tab *tab, int row)
1040 int j;
1041 int col = tab->n_col;
1042 isl_int *tr;
1043 isl_int tmp;
1045 tr = tab->mat->row[row] + 2 + tab->M;
1047 isl_int_init(tmp);
1049 for (j = tab->n_dead; j < tab->n_col; ++j) {
1050 if (tab->col_var[j] >= 0 &&
1051 (tab->col_var[j] < tab->n_param ||
1052 tab->col_var[j] >= tab->n_var - tab->n_div))
1053 continue;
1055 if (!isl_int_is_pos(tr[j]))
1056 continue;
1058 if (col == tab->n_col)
1059 col = j;
1060 else
1061 col = lexmin_col_pair(tab, row, col, j, tmp);
1062 isl_assert(tab->mat->ctx, col >= 0, goto error);
1065 isl_int_clear(tmp);
1066 return col;
1067 error:
1068 isl_int_clear(tmp);
1069 return -1;
1072 /* Return the first known violated constraint, i.e., a non-negative
1073 * contraint that currently has an either obviously negative value
1074 * or a previously determined to be negative value.
1076 * If any constraint has a negative coefficient for the big parameter,
1077 * if any, then we return one of these first.
1079 static int first_neg(struct isl_tab *tab)
1081 int row;
1083 if (tab->M)
1084 for (row = tab->n_redundant; row < tab->n_row; ++row) {
1085 if (!isl_tab_var_from_row(tab, row)->is_nonneg)
1086 continue;
1087 if (!isl_int_is_neg(tab->mat->row[row][2]))
1088 continue;
1089 if (tab->row_sign)
1090 tab->row_sign[row] = isl_tab_row_neg;
1091 return row;
1093 for (row = tab->n_redundant; row < tab->n_row; ++row) {
1094 if (!isl_tab_var_from_row(tab, row)->is_nonneg)
1095 continue;
1096 if (tab->row_sign) {
1097 if (tab->row_sign[row] == 0 &&
1098 is_obviously_neg(tab, row))
1099 tab->row_sign[row] = isl_tab_row_neg;
1100 if (tab->row_sign[row] != isl_tab_row_neg)
1101 continue;
1102 } else if (!is_obviously_neg(tab, row))
1103 continue;
1104 return row;
1106 return -1;
1109 /* Resolve all known or obviously violated constraints through pivoting.
1110 * In particular, as long as we can find any violated constraint, we
1111 * look for a pivoting column that would result in the lexicographicallly
1112 * smallest increment in the sample point. If there is no such column
1113 * then the tableau is infeasible.
1115 static struct isl_tab *restore_lexmin(struct isl_tab *tab) WARN_UNUSED;
1116 static struct isl_tab *restore_lexmin(struct isl_tab *tab)
1118 int row, col;
1120 if (!tab)
1121 return NULL;
1122 if (tab->empty)
1123 return tab;
1124 while ((row = first_neg(tab)) != -1) {
1125 col = lexmin_pivot_col(tab, row);
1126 if (col >= tab->n_col) {
1127 if (isl_tab_mark_empty(tab) < 0)
1128 goto error;
1129 return tab;
1131 if (col < 0)
1132 goto error;
1133 if (isl_tab_pivot(tab, row, col) < 0)
1134 goto error;
1136 return tab;
1137 error:
1138 isl_tab_free(tab);
1139 return NULL;
1142 /* Given a row that represents an equality, look for an appropriate
1143 * pivoting column.
1144 * In particular, if there are any non-zero coefficients among
1145 * the non-parameter variables, then we take the last of these
1146 * variables. Eliminating this variable in terms of the other
1147 * variables and/or parameters does not influence the property
1148 * that all column in the initial tableau are lexicographically
1149 * positive. The row corresponding to the eliminated variable
1150 * will only have non-zero entries below the diagonal of the
1151 * initial tableau. That is, we transform
1153 * I I
1154 * 1 into a
1155 * I I
1157 * If there is no such non-parameter variable, then we are dealing with
1158 * pure parameter equality and we pick any parameter with coefficient 1 or -1
1159 * for elimination. This will ensure that the eliminated parameter
1160 * always has an integer value whenever all the other parameters are integral.
1161 * If there is no such parameter then we return -1.
1163 static int last_var_col_or_int_par_col(struct isl_tab *tab, int row)
1165 unsigned off = 2 + tab->M;
1166 int i;
1168 for (i = tab->n_var - tab->n_div - 1; i >= 0 && i >= tab->n_param; --i) {
1169 int col;
1170 if (tab->var[i].is_row)
1171 continue;
1172 col = tab->var[i].index;
1173 if (col <= tab->n_dead)
1174 continue;
1175 if (!isl_int_is_zero(tab->mat->row[row][off + col]))
1176 return col;
1178 for (i = tab->n_dead; i < tab->n_col; ++i) {
1179 if (isl_int_is_one(tab->mat->row[row][off + i]))
1180 return i;
1181 if (isl_int_is_negone(tab->mat->row[row][off + i]))
1182 return i;
1184 return -1;
1187 /* Add an equality that is known to be valid to the tableau.
1188 * We first check if we can eliminate a variable or a parameter.
1189 * If not, we add the equality as two inequalities.
1190 * In this case, the equality was a pure parameter equality and there
1191 * is no need to resolve any constraint violations.
1193 static struct isl_tab *add_lexmin_valid_eq(struct isl_tab *tab, isl_int *eq)
1195 int i;
1196 int r;
1198 if (!tab)
1199 return NULL;
1200 r = isl_tab_add_row(tab, eq);
1201 if (r < 0)
1202 goto error;
1204 r = tab->con[r].index;
1205 i = last_var_col_or_int_par_col(tab, r);
1206 if (i < 0) {
1207 tab->con[r].is_nonneg = 1;
1208 if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r]) < 0)
1209 goto error;
1210 isl_seq_neg(eq, eq, 1 + tab->n_var);
1211 r = isl_tab_add_row(tab, eq);
1212 if (r < 0)
1213 goto error;
1214 tab->con[r].is_nonneg = 1;
1215 if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r]) < 0)
1216 goto error;
1217 } else {
1218 if (isl_tab_pivot(tab, r, i) < 0)
1219 goto error;
1220 if (isl_tab_kill_col(tab, i) < 0)
1221 goto error;
1222 tab->n_eq++;
1224 tab = restore_lexmin(tab);
1227 return tab;
1228 error:
1229 isl_tab_free(tab);
1230 return NULL;
1233 /* Check if the given row is a pure constant.
1235 static int is_constant(struct isl_tab *tab, int row)
1237 unsigned off = 2 + tab->M;
1239 return isl_seq_first_non_zero(tab->mat->row[row] + off + tab->n_dead,
1240 tab->n_col - tab->n_dead) == -1;
1243 /* Add an equality that may or may not be valid to the tableau.
1244 * If the resulting row is a pure constant, then it must be zero.
1245 * Otherwise, the resulting tableau is empty.
1247 * If the row is not a pure constant, then we add two inequalities,
1248 * each time checking that they can be satisfied.
1249 * In the end we try to use one of the two constraints to eliminate
1250 * a column.
1252 static struct isl_tab *add_lexmin_eq(struct isl_tab *tab, isl_int *eq) WARN_UNUSED;
1253 static struct isl_tab *add_lexmin_eq(struct isl_tab *tab, isl_int *eq)
1255 int r1, r2;
1256 int row;
1257 struct isl_tab_undo *snap;
1259 if (!tab)
1260 return NULL;
1261 snap = isl_tab_snap(tab);
1262 r1 = isl_tab_add_row(tab, eq);
1263 if (r1 < 0)
1264 goto error;
1265 tab->con[r1].is_nonneg = 1;
1266 if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r1]) < 0)
1267 goto error;
1269 row = tab->con[r1].index;
1270 if (is_constant(tab, row)) {
1271 if (!isl_int_is_zero(tab->mat->row[row][1]) ||
1272 (tab->M && !isl_int_is_zero(tab->mat->row[row][2]))) {
1273 if (isl_tab_mark_empty(tab) < 0)
1274 goto error;
1275 return tab;
1277 if (isl_tab_rollback(tab, snap) < 0)
1278 goto error;
1279 return tab;
1282 tab = restore_lexmin(tab);
1283 if (!tab || tab->empty)
1284 return tab;
1286 isl_seq_neg(eq, eq, 1 + tab->n_var);
1288 r2 = isl_tab_add_row(tab, eq);
1289 if (r2 < 0)
1290 goto error;
1291 tab->con[r2].is_nonneg = 1;
1292 if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r2]) < 0)
1293 goto error;
1295 tab = restore_lexmin(tab);
1296 if (!tab || tab->empty)
1297 return tab;
1299 if (!tab->con[r1].is_row) {
1300 if (isl_tab_kill_col(tab, tab->con[r1].index) < 0)
1301 goto error;
1302 } else if (!tab->con[r2].is_row) {
1303 if (isl_tab_kill_col(tab, tab->con[r2].index) < 0)
1304 goto error;
1305 } else if (isl_int_is_zero(tab->mat->row[tab->con[r1].index][1])) {
1306 unsigned off = 2 + tab->M;
1307 int i;
1308 int row = tab->con[r1].index;
1309 i = isl_seq_first_non_zero(tab->mat->row[row]+off+tab->n_dead,
1310 tab->n_col - tab->n_dead);
1311 if (i != -1) {
1312 if (isl_tab_pivot(tab, row, tab->n_dead + i) < 0)
1313 goto error;
1314 if (isl_tab_kill_col(tab, tab->n_dead + i) < 0)
1315 goto error;
1319 if (tab->bmap) {
1320 tab->bmap = isl_basic_map_add_ineq(tab->bmap, eq);
1321 if (isl_tab_push(tab, isl_tab_undo_bmap_ineq) < 0)
1322 goto error;
1323 isl_seq_neg(eq, eq, 1 + tab->n_var);
1324 tab->bmap = isl_basic_map_add_ineq(tab->bmap, eq);
1325 isl_seq_neg(eq, eq, 1 + tab->n_var);
1326 if (isl_tab_push(tab, isl_tab_undo_bmap_ineq) < 0)
1327 goto error;
1328 if (!tab->bmap)
1329 goto error;
1332 return tab;
1333 error:
1334 isl_tab_free(tab);
1335 return NULL;
1338 /* Add an inequality to the tableau, resolving violations using
1339 * restore_lexmin.
1341 static struct isl_tab *add_lexmin_ineq(struct isl_tab *tab, isl_int *ineq)
1343 int r;
1345 if (!tab)
1346 return NULL;
1347 if (tab->bmap) {
1348 tab->bmap = isl_basic_map_add_ineq(tab->bmap, ineq);
1349 if (isl_tab_push(tab, isl_tab_undo_bmap_ineq) < 0)
1350 goto error;
1351 if (!tab->bmap)
1352 goto error;
1354 r = isl_tab_add_row(tab, ineq);
1355 if (r < 0)
1356 goto error;
1357 tab->con[r].is_nonneg = 1;
1358 if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r]) < 0)
1359 goto error;
1360 if (isl_tab_row_is_redundant(tab, tab->con[r].index)) {
1361 if (isl_tab_mark_redundant(tab, tab->con[r].index) < 0)
1362 goto error;
1363 return tab;
1366 tab = restore_lexmin(tab);
1367 if (tab && !tab->empty && tab->con[r].is_row &&
1368 isl_tab_row_is_redundant(tab, tab->con[r].index))
1369 if (isl_tab_mark_redundant(tab, tab->con[r].index) < 0)
1370 goto error;
1371 return tab;
1372 error:
1373 isl_tab_free(tab);
1374 return NULL;
1377 /* Check if the coefficients of the parameters are all integral.
1379 static int integer_parameter(struct isl_tab *tab, int row)
1381 int i;
1382 int col;
1383 unsigned off = 2 + tab->M;
1385 for (i = 0; i < tab->n_param; ++i) {
1386 /* Eliminated parameter */
1387 if (tab->var[i].is_row)
1388 continue;
1389 col = tab->var[i].index;
1390 if (!isl_int_is_divisible_by(tab->mat->row[row][off + col],
1391 tab->mat->row[row][0]))
1392 return 0;
1394 for (i = 0; i < tab->n_div; ++i) {
1395 if (tab->var[tab->n_var - tab->n_div + i].is_row)
1396 continue;
1397 col = tab->var[tab->n_var - tab->n_div + i].index;
1398 if (!isl_int_is_divisible_by(tab->mat->row[row][off + col],
1399 tab->mat->row[row][0]))
1400 return 0;
1402 return 1;
1405 /* Check if the coefficients of the non-parameter variables are all integral.
1407 static int integer_variable(struct isl_tab *tab, int row)
1409 int i;
1410 unsigned off = 2 + tab->M;
1412 for (i = tab->n_dead; i < tab->n_col; ++i) {
1413 if (tab->col_var[i] >= 0 &&
1414 (tab->col_var[i] < tab->n_param ||
1415 tab->col_var[i] >= tab->n_var - tab->n_div))
1416 continue;
1417 if (!isl_int_is_divisible_by(tab->mat->row[row][off + i],
1418 tab->mat->row[row][0]))
1419 return 0;
1421 return 1;
1424 /* Check if the constant term is integral.
1426 static int integer_constant(struct isl_tab *tab, int row)
1428 return isl_int_is_divisible_by(tab->mat->row[row][1],
1429 tab->mat->row[row][0]);
1432 #define I_CST 1 << 0
1433 #define I_PAR 1 << 1
1434 #define I_VAR 1 << 2
1436 /* Check for next (non-parameter) variable after "var" (first if var == -1)
1437 * that is non-integer and therefore requires a cut and return
1438 * the index of the variable.
1439 * For parametric tableaus, there are three parts in a row,
1440 * the constant, the coefficients of the parameters and the rest.
1441 * For each part, we check whether the coefficients in that part
1442 * are all integral and if so, set the corresponding flag in *f.
1443 * If the constant and the parameter part are integral, then the
1444 * current sample value is integral and no cut is required
1445 * (irrespective of whether the variable part is integral).
1447 static int next_non_integer_var(struct isl_tab *tab, int var, int *f)
1449 var = var < 0 ? tab->n_param : var + 1;
1451 for (; var < tab->n_var - tab->n_div; ++var) {
1452 int flags = 0;
1453 int row;
1454 if (!tab->var[var].is_row)
1455 continue;
1456 row = tab->var[var].index;
1457 if (integer_constant(tab, row))
1458 ISL_FL_SET(flags, I_CST);
1459 if (integer_parameter(tab, row))
1460 ISL_FL_SET(flags, I_PAR);
1461 if (ISL_FL_ISSET(flags, I_CST) && ISL_FL_ISSET(flags, I_PAR))
1462 continue;
1463 if (integer_variable(tab, row))
1464 ISL_FL_SET(flags, I_VAR);
1465 *f = flags;
1466 return var;
1468 return -1;
1471 /* Check for first (non-parameter) variable that is non-integer and
1472 * therefore requires a cut and return the corresponding row.
1473 * For parametric tableaus, there are three parts in a row,
1474 * the constant, the coefficients of the parameters and the rest.
1475 * For each part, we check whether the coefficients in that part
1476 * are all integral and if so, set the corresponding flag in *f.
1477 * If the constant and the parameter part are integral, then the
1478 * current sample value is integral and no cut is required
1479 * (irrespective of whether the variable part is integral).
1481 static int first_non_integer_row(struct isl_tab *tab, int *f)
1483 int var = next_non_integer_var(tab, -1, f);
1485 return var < 0 ? -1 : tab->var[var].index;
1488 /* Add a (non-parametric) cut to cut away the non-integral sample
1489 * value of the given row.
1491 * If the row is given by
1493 * m r = f + \sum_i a_i y_i
1495 * then the cut is
1497 * c = - {-f/m} + \sum_i {a_i/m} y_i >= 0
1499 * The big parameter, if any, is ignored, since it is assumed to be big
1500 * enough to be divisible by any integer.
1501 * If the tableau is actually a parametric tableau, then this function
1502 * is only called when all coefficients of the parameters are integral.
1503 * The cut therefore has zero coefficients for the parameters.
1505 * The current value is known to be negative, so row_sign, if it
1506 * exists, is set accordingly.
1508 * Return the row of the cut or -1.
1510 static int add_cut(struct isl_tab *tab, int row)
1512 int i;
1513 int r;
1514 isl_int *r_row;
1515 unsigned off = 2 + tab->M;
1517 if (isl_tab_extend_cons(tab, 1) < 0)
1518 return -1;
1519 r = isl_tab_allocate_con(tab);
1520 if (r < 0)
1521 return -1;
1523 r_row = tab->mat->row[tab->con[r].index];
1524 isl_int_set(r_row[0], tab->mat->row[row][0]);
1525 isl_int_neg(r_row[1], tab->mat->row[row][1]);
1526 isl_int_fdiv_r(r_row[1], r_row[1], tab->mat->row[row][0]);
1527 isl_int_neg(r_row[1], r_row[1]);
1528 if (tab->M)
1529 isl_int_set_si(r_row[2], 0);
1530 for (i = 0; i < tab->n_col; ++i)
1531 isl_int_fdiv_r(r_row[off + i],
1532 tab->mat->row[row][off + i], tab->mat->row[row][0]);
1534 tab->con[r].is_nonneg = 1;
1535 if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r]) < 0)
1536 return -1;
1537 if (tab->row_sign)
1538 tab->row_sign[tab->con[r].index] = isl_tab_row_neg;
1540 return tab->con[r].index;
1543 /* Given a non-parametric tableau, add cuts until an integer
1544 * sample point is obtained or until the tableau is determined
1545 * to be integer infeasible.
1546 * As long as there is any non-integer value in the sample point,
1547 * we add appropriate cuts, if possible, for each of these
1548 * non-integer values and then resolve the violated
1549 * cut constraints using restore_lexmin.
1550 * If one of the corresponding rows is equal to an integral
1551 * combination of variables/constraints plus a non-integral constant,
1552 * then there is no way to obtain an integer point and we return
1553 * a tableau that is marked empty.
1555 static struct isl_tab *cut_to_integer_lexmin(struct isl_tab *tab)
1557 int var;
1558 int row;
1559 int flags;
1561 if (!tab)
1562 return NULL;
1563 if (tab->empty)
1564 return tab;
1566 while ((var = next_non_integer_var(tab, -1, &flags)) != -1) {
1567 do {
1568 if (ISL_FL_ISSET(flags, I_VAR)) {
1569 if (isl_tab_mark_empty(tab) < 0)
1570 goto error;
1571 return tab;
1573 row = tab->var[var].index;
1574 row = add_cut(tab, row);
1575 if (row < 0)
1576 goto error;
1577 } while ((var = next_non_integer_var(tab, var, &flags)) != -1);
1578 tab = restore_lexmin(tab);
1579 if (!tab || tab->empty)
1580 break;
1582 return tab;
1583 error:
1584 isl_tab_free(tab);
1585 return NULL;
1588 /* Check whether all the currently active samples also satisfy the inequality
1589 * "ineq" (treated as an equality if eq is set).
1590 * Remove those samples that do not.
1592 static struct isl_tab *check_samples(struct isl_tab *tab, isl_int *ineq, int eq)
1594 int i;
1595 isl_int v;
1597 if (!tab)
1598 return NULL;
1600 isl_assert(tab->mat->ctx, tab->bmap, goto error);
1601 isl_assert(tab->mat->ctx, tab->samples, goto error);
1602 isl_assert(tab->mat->ctx, tab->samples->n_col == 1 + tab->n_var, goto error);
1604 isl_int_init(v);
1605 for (i = tab->n_outside; i < tab->n_sample; ++i) {
1606 int sgn;
1607 isl_seq_inner_product(ineq, tab->samples->row[i],
1608 1 + tab->n_var, &v);
1609 sgn = isl_int_sgn(v);
1610 if (eq ? (sgn == 0) : (sgn >= 0))
1611 continue;
1612 tab = isl_tab_drop_sample(tab, i);
1613 if (!tab)
1614 break;
1616 isl_int_clear(v);
1618 return tab;
1619 error:
1620 isl_tab_free(tab);
1621 return NULL;
1624 /* Check whether the sample value of the tableau is finite,
1625 * i.e., either the tableau does not use a big parameter, or
1626 * all values of the variables are equal to the big parameter plus
1627 * some constant. This constant is the actual sample value.
1629 static int sample_is_finite(struct isl_tab *tab)
1631 int i;
1633 if (!tab->M)
1634 return 1;
1636 for (i = 0; i < tab->n_var; ++i) {
1637 int row;
1638 if (!tab->var[i].is_row)
1639 return 0;
1640 row = tab->var[i].index;
1641 if (isl_int_ne(tab->mat->row[row][0], tab->mat->row[row][2]))
1642 return 0;
1644 return 1;
1647 /* Check if the context tableau of sol has any integer points.
1648 * Leave tab in empty state if no integer point can be found.
1649 * If an integer point can be found and if moreover it is finite,
1650 * then it is added to the list of sample values.
1652 * This function is only called when none of the currently active sample
1653 * values satisfies the most recently added constraint.
1655 static struct isl_tab *check_integer_feasible(struct isl_tab *tab)
1657 struct isl_tab_undo *snap;
1658 int feasible;
1660 if (!tab)
1661 return NULL;
1663 snap = isl_tab_snap(tab);
1664 if (isl_tab_push_basis(tab) < 0)
1665 goto error;
1667 tab = cut_to_integer_lexmin(tab);
1668 if (!tab)
1669 goto error;
1671 if (!tab->empty && sample_is_finite(tab)) {
1672 struct isl_vec *sample;
1674 sample = isl_tab_get_sample_value(tab);
1676 tab = isl_tab_add_sample(tab, sample);
1679 if (!tab->empty && isl_tab_rollback(tab, snap) < 0)
1680 goto error;
1682 return tab;
1683 error:
1684 isl_tab_free(tab);
1685 return NULL;
1688 /* Check if any of the currently active sample values satisfies
1689 * the inequality "ineq" (an equality if eq is set).
1691 static int tab_has_valid_sample(struct isl_tab *tab, isl_int *ineq, int eq)
1693 int i;
1694 isl_int v;
1696 if (!tab)
1697 return -1;
1699 isl_assert(tab->mat->ctx, tab->bmap, return -1);
1700 isl_assert(tab->mat->ctx, tab->samples, return -1);
1701 isl_assert(tab->mat->ctx, tab->samples->n_col == 1 + tab->n_var, return -1);
1703 isl_int_init(v);
1704 for (i = tab->n_outside; i < tab->n_sample; ++i) {
1705 int sgn;
1706 isl_seq_inner_product(ineq, tab->samples->row[i],
1707 1 + tab->n_var, &v);
1708 sgn = isl_int_sgn(v);
1709 if (eq ? (sgn == 0) : (sgn >= 0))
1710 break;
1712 isl_int_clear(v);
1714 return i < tab->n_sample;
1717 /* Add a div specifed by "div" to the tableau "tab" and return
1718 * 1 if the div is obviously non-negative.
1720 static int context_tab_add_div(struct isl_tab *tab, struct isl_vec *div,
1721 int (*add_ineq)(void *user, isl_int *), void *user)
1723 int i;
1724 int r;
1725 struct isl_mat *samples;
1726 int nonneg;
1728 r = isl_tab_add_div(tab, div, add_ineq, user);
1729 if (r < 0)
1730 return -1;
1731 nonneg = tab->var[r].is_nonneg;
1732 tab->var[r].frozen = 1;
1734 samples = isl_mat_extend(tab->samples,
1735 tab->n_sample, 1 + tab->n_var);
1736 tab->samples = samples;
1737 if (!samples)
1738 return -1;
1739 for (i = tab->n_outside; i < samples->n_row; ++i) {
1740 isl_seq_inner_product(div->el + 1, samples->row[i],
1741 div->size - 1, &samples->row[i][samples->n_col - 1]);
1742 isl_int_fdiv_q(samples->row[i][samples->n_col - 1],
1743 samples->row[i][samples->n_col - 1], div->el[0]);
1746 return nonneg;
1749 /* Add a div specified by "div" to both the main tableau and
1750 * the context tableau. In case of the main tableau, we only
1751 * need to add an extra div. In the context tableau, we also
1752 * need to express the meaning of the div.
1753 * Return the index of the div or -1 if anything went wrong.
1755 static int add_div(struct isl_tab *tab, struct isl_context *context,
1756 struct isl_vec *div)
1758 int r;
1759 int nonneg;
1761 if ((nonneg = context->op->add_div(context, div)) < 0)
1762 goto error;
1764 if (!context->op->is_ok(context))
1765 goto error;
1767 if (isl_tab_extend_vars(tab, 1) < 0)
1768 goto error;
1769 r = isl_tab_allocate_var(tab);
1770 if (r < 0)
1771 goto error;
1772 if (nonneg)
1773 tab->var[r].is_nonneg = 1;
1774 tab->var[r].frozen = 1;
1775 tab->n_div++;
1777 return tab->n_div - 1;
1778 error:
1779 context->op->invalidate(context);
1780 return -1;
1783 static int find_div(struct isl_tab *tab, isl_int *div, isl_int denom)
1785 int i;
1786 unsigned total = isl_basic_map_total_dim(tab->bmap);
1788 for (i = 0; i < tab->bmap->n_div; ++i) {
1789 if (isl_int_ne(tab->bmap->div[i][0], denom))
1790 continue;
1791 if (!isl_seq_eq(tab->bmap->div[i] + 1, div, total))
1792 continue;
1793 return i;
1795 return -1;
1798 /* Return the index of a div that corresponds to "div".
1799 * We first check if we already have such a div and if not, we create one.
1801 static int get_div(struct isl_tab *tab, struct isl_context *context,
1802 struct isl_vec *div)
1804 int d;
1805 struct isl_tab *context_tab = context->op->peek_tab(context);
1807 if (!context_tab)
1808 return -1;
1810 d = find_div(context_tab, div->el + 1, div->el[0]);
1811 if (d != -1)
1812 return d;
1814 return add_div(tab, context, div);
1817 /* Add a parametric cut to cut away the non-integral sample value
1818 * of the give row.
1819 * Let a_i be the coefficients of the constant term and the parameters
1820 * and let b_i be the coefficients of the variables or constraints
1821 * in basis of the tableau.
1822 * Let q be the div q = floor(\sum_i {-a_i} y_i).
1824 * The cut is expressed as
1826 * c = \sum_i -{-a_i} y_i + \sum_i {b_i} x_i + q >= 0
1828 * If q did not already exist in the context tableau, then it is added first.
1829 * If q is in a column of the main tableau then the "+ q" can be accomplished
1830 * by setting the corresponding entry to the denominator of the constraint.
1831 * If q happens to be in a row of the main tableau, then the corresponding
1832 * row needs to be added instead (taking care of the denominators).
1833 * Note that this is very unlikely, but perhaps not entirely impossible.
1835 * The current value of the cut is known to be negative (or at least
1836 * non-positive), so row_sign is set accordingly.
1838 * Return the row of the cut or -1.
1840 static int add_parametric_cut(struct isl_tab *tab, int row,
1841 struct isl_context *context)
1843 struct isl_vec *div;
1844 int d;
1845 int i;
1846 int r;
1847 isl_int *r_row;
1848 int col;
1849 int n;
1850 unsigned off = 2 + tab->M;
1852 if (!context)
1853 return -1;
1855 div = get_row_parameter_div(tab, row);
1856 if (!div)
1857 return -1;
1859 n = tab->n_div;
1860 d = context->op->get_div(context, tab, div);
1861 if (d < 0)
1862 return -1;
1864 if (isl_tab_extend_cons(tab, 1) < 0)
1865 return -1;
1866 r = isl_tab_allocate_con(tab);
1867 if (r < 0)
1868 return -1;
1870 r_row = tab->mat->row[tab->con[r].index];
1871 isl_int_set(r_row[0], tab->mat->row[row][0]);
1872 isl_int_neg(r_row[1], tab->mat->row[row][1]);
1873 isl_int_fdiv_r(r_row[1], r_row[1], tab->mat->row[row][0]);
1874 isl_int_neg(r_row[1], r_row[1]);
1875 if (tab->M)
1876 isl_int_set_si(r_row[2], 0);
1877 for (i = 0; i < tab->n_param; ++i) {
1878 if (tab->var[i].is_row)
1879 continue;
1880 col = tab->var[i].index;
1881 isl_int_neg(r_row[off + col], tab->mat->row[row][off + col]);
1882 isl_int_fdiv_r(r_row[off + col], r_row[off + col],
1883 tab->mat->row[row][0]);
1884 isl_int_neg(r_row[off + col], r_row[off + col]);
1886 for (i = 0; i < tab->n_div; ++i) {
1887 if (tab->var[tab->n_var - tab->n_div + i].is_row)
1888 continue;
1889 col = tab->var[tab->n_var - tab->n_div + i].index;
1890 isl_int_neg(r_row[off + col], tab->mat->row[row][off + col]);
1891 isl_int_fdiv_r(r_row[off + col], r_row[off + col],
1892 tab->mat->row[row][0]);
1893 isl_int_neg(r_row[off + col], r_row[off + col]);
1895 for (i = 0; i < tab->n_col; ++i) {
1896 if (tab->col_var[i] >= 0 &&
1897 (tab->col_var[i] < tab->n_param ||
1898 tab->col_var[i] >= tab->n_var - tab->n_div))
1899 continue;
1900 isl_int_fdiv_r(r_row[off + i],
1901 tab->mat->row[row][off + i], tab->mat->row[row][0]);
1903 if (tab->var[tab->n_var - tab->n_div + d].is_row) {
1904 isl_int gcd;
1905 int d_row = tab->var[tab->n_var - tab->n_div + d].index;
1906 isl_int_init(gcd);
1907 isl_int_gcd(gcd, tab->mat->row[d_row][0], r_row[0]);
1908 isl_int_divexact(r_row[0], r_row[0], gcd);
1909 isl_int_divexact(gcd, tab->mat->row[d_row][0], gcd);
1910 isl_seq_combine(r_row + 1, gcd, r_row + 1,
1911 r_row[0], tab->mat->row[d_row] + 1,
1912 off - 1 + tab->n_col);
1913 isl_int_mul(r_row[0], r_row[0], tab->mat->row[d_row][0]);
1914 isl_int_clear(gcd);
1915 } else {
1916 col = tab->var[tab->n_var - tab->n_div + d].index;
1917 isl_int_set(r_row[off + col], tab->mat->row[row][0]);
1920 tab->con[r].is_nonneg = 1;
1921 if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r]) < 0)
1922 return -1;
1923 if (tab->row_sign)
1924 tab->row_sign[tab->con[r].index] = isl_tab_row_neg;
1926 isl_vec_free(div);
1928 row = tab->con[r].index;
1930 if (d >= n && context->op->detect_equalities(context, tab) < 0)
1931 return -1;
1933 return row;
1936 /* Construct a tableau for bmap that can be used for computing
1937 * the lexicographic minimum (or maximum) of bmap.
1938 * If not NULL, then dom is the domain where the minimum
1939 * should be computed. In this case, we set up a parametric
1940 * tableau with row signs (initialized to "unknown").
1941 * If M is set, then the tableau will use a big parameter.
1942 * If max is set, then a maximum should be computed instead of a minimum.
1943 * This means that for each variable x, the tableau will contain the variable
1944 * x' = M - x, rather than x' = M + x. This in turn means that the coefficient
1945 * of the variables in all constraints are negated prior to adding them
1946 * to the tableau.
1948 static struct isl_tab *tab_for_lexmin(struct isl_basic_map *bmap,
1949 struct isl_basic_set *dom, unsigned M, int max)
1951 int i;
1952 struct isl_tab *tab;
1954 tab = isl_tab_alloc(bmap->ctx, 2 * bmap->n_eq + bmap->n_ineq + 1,
1955 isl_basic_map_total_dim(bmap), M);
1956 if (!tab)
1957 return NULL;
1959 tab->rational = ISL_F_ISSET(bmap, ISL_BASIC_MAP_RATIONAL);
1960 if (dom) {
1961 tab->n_param = isl_basic_set_total_dim(dom) - dom->n_div;
1962 tab->n_div = dom->n_div;
1963 tab->row_sign = isl_calloc_array(bmap->ctx,
1964 enum isl_tab_row_sign, tab->mat->n_row);
1965 if (!tab->row_sign)
1966 goto error;
1968 if (ISL_F_ISSET(bmap, ISL_BASIC_MAP_EMPTY)) {
1969 if (isl_tab_mark_empty(tab) < 0)
1970 goto error;
1971 return tab;
1974 for (i = tab->n_param; i < tab->n_var - tab->n_div; ++i) {
1975 tab->var[i].is_nonneg = 1;
1976 tab->var[i].frozen = 1;
1978 for (i = 0; i < bmap->n_eq; ++i) {
1979 if (max)
1980 isl_seq_neg(bmap->eq[i] + 1 + tab->n_param,
1981 bmap->eq[i] + 1 + tab->n_param,
1982 tab->n_var - tab->n_param - tab->n_div);
1983 tab = add_lexmin_valid_eq(tab, bmap->eq[i]);
1984 if (max)
1985 isl_seq_neg(bmap->eq[i] + 1 + tab->n_param,
1986 bmap->eq[i] + 1 + tab->n_param,
1987 tab->n_var - tab->n_param - tab->n_div);
1988 if (!tab || tab->empty)
1989 return tab;
1991 for (i = 0; i < bmap->n_ineq; ++i) {
1992 if (max)
1993 isl_seq_neg(bmap->ineq[i] + 1 + tab->n_param,
1994 bmap->ineq[i] + 1 + tab->n_param,
1995 tab->n_var - tab->n_param - tab->n_div);
1996 tab = add_lexmin_ineq(tab, bmap->ineq[i]);
1997 if (max)
1998 isl_seq_neg(bmap->ineq[i] + 1 + tab->n_param,
1999 bmap->ineq[i] + 1 + tab->n_param,
2000 tab->n_var - tab->n_param - tab->n_div);
2001 if (!tab || tab->empty)
2002 return tab;
2004 return tab;
2005 error:
2006 isl_tab_free(tab);
2007 return NULL;
2010 /* Given a main tableau where more than one row requires a split,
2011 * determine and return the "best" row to split on.
2013 * Given two rows in the main tableau, if the inequality corresponding
2014 * to the first row is redundant with respect to that of the second row
2015 * in the current tableau, then it is better to split on the second row,
2016 * since in the positive part, both row will be positive.
2017 * (In the negative part a pivot will have to be performed and just about
2018 * anything can happen to the sign of the other row.)
2020 * As a simple heuristic, we therefore select the row that makes the most
2021 * of the other rows redundant.
2023 * Perhaps it would also be useful to look at the number of constraints
2024 * that conflict with any given constraint.
2026 static int best_split(struct isl_tab *tab, struct isl_tab *context_tab)
2028 struct isl_tab_undo *snap;
2029 int split;
2030 int row;
2031 int best = -1;
2032 int best_r;
2034 if (isl_tab_extend_cons(context_tab, 2) < 0)
2035 return -1;
2037 snap = isl_tab_snap(context_tab);
2039 for (split = tab->n_redundant; split < tab->n_row; ++split) {
2040 struct isl_tab_undo *snap2;
2041 struct isl_vec *ineq = NULL;
2042 int r = 0;
2043 int ok;
2045 if (!isl_tab_var_from_row(tab, split)->is_nonneg)
2046 continue;
2047 if (tab->row_sign[split] != isl_tab_row_any)
2048 continue;
2050 ineq = get_row_parameter_ineq(tab, split);
2051 if (!ineq)
2052 return -1;
2053 ok = isl_tab_add_ineq(context_tab, ineq->el) >= 0;
2054 isl_vec_free(ineq);
2055 if (!ok)
2056 return -1;
2058 snap2 = isl_tab_snap(context_tab);
2060 for (row = tab->n_redundant; row < tab->n_row; ++row) {
2061 struct isl_tab_var *var;
2063 if (row == split)
2064 continue;
2065 if (!isl_tab_var_from_row(tab, row)->is_nonneg)
2066 continue;
2067 if (tab->row_sign[row] != isl_tab_row_any)
2068 continue;
2070 ineq = get_row_parameter_ineq(tab, row);
2071 if (!ineq)
2072 return -1;
2073 ok = isl_tab_add_ineq(context_tab, ineq->el) >= 0;
2074 isl_vec_free(ineq);
2075 if (!ok)
2076 return -1;
2077 var = &context_tab->con[context_tab->n_con - 1];
2078 if (!context_tab->empty &&
2079 !isl_tab_min_at_most_neg_one(context_tab, var))
2080 r++;
2081 if (isl_tab_rollback(context_tab, snap2) < 0)
2082 return -1;
2084 if (best == -1 || r > best_r) {
2085 best = split;
2086 best_r = r;
2088 if (isl_tab_rollback(context_tab, snap) < 0)
2089 return -1;
2092 return best;
2095 static struct isl_basic_set *context_lex_peek_basic_set(
2096 struct isl_context *context)
2098 struct isl_context_lex *clex = (struct isl_context_lex *)context;
2099 if (!clex->tab)
2100 return NULL;
2101 return isl_tab_peek_bset(clex->tab);
2104 static struct isl_tab *context_lex_peek_tab(struct isl_context *context)
2106 struct isl_context_lex *clex = (struct isl_context_lex *)context;
2107 return clex->tab;
2110 static void context_lex_extend(struct isl_context *context, int n)
2112 struct isl_context_lex *clex = (struct isl_context_lex *)context;
2113 if (!clex->tab)
2114 return;
2115 if (isl_tab_extend_cons(clex->tab, n) >= 0)
2116 return;
2117 isl_tab_free(clex->tab);
2118 clex->tab = NULL;
2121 static void context_lex_add_eq(struct isl_context *context, isl_int *eq,
2122 int check, int update)
2124 struct isl_context_lex *clex = (struct isl_context_lex *)context;
2125 if (isl_tab_extend_cons(clex->tab, 2) < 0)
2126 goto error;
2127 clex->tab = add_lexmin_eq(clex->tab, eq);
2128 if (check) {
2129 int v = tab_has_valid_sample(clex->tab, eq, 1);
2130 if (v < 0)
2131 goto error;
2132 if (!v)
2133 clex->tab = check_integer_feasible(clex->tab);
2135 if (update)
2136 clex->tab = check_samples(clex->tab, eq, 1);
2137 return;
2138 error:
2139 isl_tab_free(clex->tab);
2140 clex->tab = NULL;
2143 static void context_lex_add_ineq(struct isl_context *context, isl_int *ineq,
2144 int check, int update)
2146 struct isl_context_lex *clex = (struct isl_context_lex *)context;
2147 if (isl_tab_extend_cons(clex->tab, 1) < 0)
2148 goto error;
2149 clex->tab = add_lexmin_ineq(clex->tab, ineq);
2150 if (check) {
2151 int v = tab_has_valid_sample(clex->tab, ineq, 0);
2152 if (v < 0)
2153 goto error;
2154 if (!v)
2155 clex->tab = check_integer_feasible(clex->tab);
2157 if (update)
2158 clex->tab = check_samples(clex->tab, ineq, 0);
2159 return;
2160 error:
2161 isl_tab_free(clex->tab);
2162 clex->tab = NULL;
2165 static int context_lex_add_ineq_wrap(void *user, isl_int *ineq)
2167 struct isl_context *context = (struct isl_context *)user;
2168 context_lex_add_ineq(context, ineq, 0, 0);
2169 return context->op->is_ok(context) ? 0 : -1;
2172 /* Check which signs can be obtained by "ineq" on all the currently
2173 * active sample values. See row_sign for more information.
2175 static enum isl_tab_row_sign tab_ineq_sign(struct isl_tab *tab, isl_int *ineq,
2176 int strict)
2178 int i;
2179 int sgn;
2180 isl_int tmp;
2181 int res = isl_tab_row_unknown;
2183 isl_assert(tab->mat->ctx, tab->samples, return 0);
2184 isl_assert(tab->mat->ctx, tab->samples->n_col == 1 + tab->n_var, return 0);
2186 isl_int_init(tmp);
2187 for (i = tab->n_outside; i < tab->n_sample; ++i) {
2188 isl_seq_inner_product(tab->samples->row[i], ineq,
2189 1 + tab->n_var, &tmp);
2190 sgn = isl_int_sgn(tmp);
2191 if (sgn > 0 || (sgn == 0 && strict)) {
2192 if (res == isl_tab_row_unknown)
2193 res = isl_tab_row_pos;
2194 if (res == isl_tab_row_neg)
2195 res = isl_tab_row_any;
2197 if (sgn < 0) {
2198 if (res == isl_tab_row_unknown)
2199 res = isl_tab_row_neg;
2200 if (res == isl_tab_row_pos)
2201 res = isl_tab_row_any;
2203 if (res == isl_tab_row_any)
2204 break;
2206 isl_int_clear(tmp);
2208 return res;
2211 static enum isl_tab_row_sign context_lex_ineq_sign(struct isl_context *context,
2212 isl_int *ineq, int strict)
2214 struct isl_context_lex *clex = (struct isl_context_lex *)context;
2215 return tab_ineq_sign(clex->tab, ineq, strict);
2218 /* Check whether "ineq" can be added to the tableau without rendering
2219 * it infeasible.
2221 static int context_lex_test_ineq(struct isl_context *context, isl_int *ineq)
2223 struct isl_context_lex *clex = (struct isl_context_lex *)context;
2224 struct isl_tab_undo *snap;
2225 int feasible;
2227 if (!clex->tab)
2228 return -1;
2230 if (isl_tab_extend_cons(clex->tab, 1) < 0)
2231 return -1;
2233 snap = isl_tab_snap(clex->tab);
2234 if (isl_tab_push_basis(clex->tab) < 0)
2235 return -1;
2236 clex->tab = add_lexmin_ineq(clex->tab, ineq);
2237 clex->tab = check_integer_feasible(clex->tab);
2238 if (!clex->tab)
2239 return -1;
2240 feasible = !clex->tab->empty;
2241 if (isl_tab_rollback(clex->tab, snap) < 0)
2242 return -1;
2244 return feasible;
2247 static int context_lex_get_div(struct isl_context *context, struct isl_tab *tab,
2248 struct isl_vec *div)
2250 return get_div(tab, context, div);
2253 static int context_lex_add_div(struct isl_context *context, struct isl_vec *div)
2255 struct isl_context_lex *clex = (struct isl_context_lex *)context;
2256 return context_tab_add_div(clex->tab, div,
2257 context_lex_add_ineq_wrap, context);
2260 static int context_lex_detect_equalities(struct isl_context *context,
2261 struct isl_tab *tab)
2263 return 0;
2266 static int context_lex_best_split(struct isl_context *context,
2267 struct isl_tab *tab)
2269 struct isl_context_lex *clex = (struct isl_context_lex *)context;
2270 struct isl_tab_undo *snap;
2271 int r;
2273 snap = isl_tab_snap(clex->tab);
2274 if (isl_tab_push_basis(clex->tab) < 0)
2275 return -1;
2276 r = best_split(tab, clex->tab);
2278 if (isl_tab_rollback(clex->tab, snap) < 0)
2279 return -1;
2281 return r;
2284 static int context_lex_is_empty(struct isl_context *context)
2286 struct isl_context_lex *clex = (struct isl_context_lex *)context;
2287 if (!clex->tab)
2288 return -1;
2289 return clex->tab->empty;
2292 static void *context_lex_save(struct isl_context *context)
2294 struct isl_context_lex *clex = (struct isl_context_lex *)context;
2295 struct isl_tab_undo *snap;
2297 snap = isl_tab_snap(clex->tab);
2298 if (isl_tab_push_basis(clex->tab) < 0)
2299 return NULL;
2300 if (isl_tab_save_samples(clex->tab) < 0)
2301 return NULL;
2303 return snap;
2306 static void context_lex_restore(struct isl_context *context, void *save)
2308 struct isl_context_lex *clex = (struct isl_context_lex *)context;
2309 if (isl_tab_rollback(clex->tab, (struct isl_tab_undo *)save) < 0) {
2310 isl_tab_free(clex->tab);
2311 clex->tab = NULL;
2315 static int context_lex_is_ok(struct isl_context *context)
2317 struct isl_context_lex *clex = (struct isl_context_lex *)context;
2318 return !!clex->tab;
2321 /* For each variable in the context tableau, check if the variable can
2322 * only attain non-negative values. If so, mark the parameter as non-negative
2323 * in the main tableau. This allows for a more direct identification of some
2324 * cases of violated constraints.
2326 static struct isl_tab *tab_detect_nonnegative_parameters(struct isl_tab *tab,
2327 struct isl_tab *context_tab)
2329 int i;
2330 struct isl_tab_undo *snap;
2331 struct isl_vec *ineq = NULL;
2332 struct isl_tab_var *var;
2333 int n;
2335 if (context_tab->n_var == 0)
2336 return tab;
2338 ineq = isl_vec_alloc(tab->mat->ctx, 1 + context_tab->n_var);
2339 if (!ineq)
2340 goto error;
2342 if (isl_tab_extend_cons(context_tab, 1) < 0)
2343 goto error;
2345 snap = isl_tab_snap(context_tab);
2347 n = 0;
2348 isl_seq_clr(ineq->el, ineq->size);
2349 for (i = 0; i < context_tab->n_var; ++i) {
2350 isl_int_set_si(ineq->el[1 + i], 1);
2351 if (isl_tab_add_ineq(context_tab, ineq->el) < 0)
2352 goto error;
2353 var = &context_tab->con[context_tab->n_con - 1];
2354 if (!context_tab->empty &&
2355 !isl_tab_min_at_most_neg_one(context_tab, var)) {
2356 int j = i;
2357 if (i >= tab->n_param)
2358 j = i - tab->n_param + tab->n_var - tab->n_div;
2359 tab->var[j].is_nonneg = 1;
2360 n++;
2362 isl_int_set_si(ineq->el[1 + i], 0);
2363 if (isl_tab_rollback(context_tab, snap) < 0)
2364 goto error;
2367 if (context_tab->M && n == context_tab->n_var) {
2368 context_tab->mat = isl_mat_drop_cols(context_tab->mat, 2, 1);
2369 context_tab->M = 0;
2372 isl_vec_free(ineq);
2373 return tab;
2374 error:
2375 isl_vec_free(ineq);
2376 isl_tab_free(tab);
2377 return NULL;
2380 static struct isl_tab *context_lex_detect_nonnegative_parameters(
2381 struct isl_context *context, struct isl_tab *tab)
2383 struct isl_context_lex *clex = (struct isl_context_lex *)context;
2384 struct isl_tab_undo *snap;
2386 snap = isl_tab_snap(clex->tab);
2387 if (isl_tab_push_basis(clex->tab) < 0)
2388 goto error;
2390 tab = tab_detect_nonnegative_parameters(tab, clex->tab);
2392 if (isl_tab_rollback(clex->tab, snap) < 0)
2393 goto error;
2395 return tab;
2396 error:
2397 isl_tab_free(tab);
2398 return NULL;
2401 static void context_lex_invalidate(struct isl_context *context)
2403 struct isl_context_lex *clex = (struct isl_context_lex *)context;
2404 isl_tab_free(clex->tab);
2405 clex->tab = NULL;
2408 static void context_lex_free(struct isl_context *context)
2410 struct isl_context_lex *clex = (struct isl_context_lex *)context;
2411 isl_tab_free(clex->tab);
2412 free(clex);
2415 struct isl_context_op isl_context_lex_op = {
2416 context_lex_detect_nonnegative_parameters,
2417 context_lex_peek_basic_set,
2418 context_lex_peek_tab,
2419 context_lex_add_eq,
2420 context_lex_add_ineq,
2421 context_lex_ineq_sign,
2422 context_lex_test_ineq,
2423 context_lex_get_div,
2424 context_lex_add_div,
2425 context_lex_detect_equalities,
2426 context_lex_best_split,
2427 context_lex_is_empty,
2428 context_lex_is_ok,
2429 context_lex_save,
2430 context_lex_restore,
2431 context_lex_invalidate,
2432 context_lex_free,
2435 static struct isl_tab *context_tab_for_lexmin(struct isl_basic_set *bset)
2437 struct isl_tab *tab;
2439 bset = isl_basic_set_cow(bset);
2440 if (!bset)
2441 return NULL;
2442 tab = tab_for_lexmin((struct isl_basic_map *)bset, NULL, 1, 0);
2443 if (!tab)
2444 goto error;
2445 if (isl_tab_track_bset(tab, bset) < 0)
2446 goto error;
2447 tab = isl_tab_init_samples(tab);
2448 return tab;
2449 error:
2450 isl_basic_set_free(bset);
2451 return NULL;
2454 static struct isl_context *isl_context_lex_alloc(struct isl_basic_set *dom)
2456 struct isl_context_lex *clex;
2458 if (!dom)
2459 return NULL;
2461 clex = isl_alloc_type(dom->ctx, struct isl_context_lex);
2462 if (!clex)
2463 return NULL;
2465 clex->context.op = &isl_context_lex_op;
2467 clex->tab = context_tab_for_lexmin(isl_basic_set_copy(dom));
2468 clex->tab = restore_lexmin(clex->tab);
2469 clex->tab = check_integer_feasible(clex->tab);
2470 if (!clex->tab)
2471 goto error;
2473 return &clex->context;
2474 error:
2475 clex->context.op->free(&clex->context);
2476 return NULL;
2479 struct isl_context_gbr {
2480 struct isl_context context;
2481 struct isl_tab *tab;
2482 struct isl_tab *shifted;
2483 struct isl_tab *cone;
2486 static struct isl_tab *context_gbr_detect_nonnegative_parameters(
2487 struct isl_context *context, struct isl_tab *tab)
2489 struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
2490 return tab_detect_nonnegative_parameters(tab, cgbr->tab);
2493 static struct isl_basic_set *context_gbr_peek_basic_set(
2494 struct isl_context *context)
2496 struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
2497 if (!cgbr->tab)
2498 return NULL;
2499 return isl_tab_peek_bset(cgbr->tab);
2502 static struct isl_tab *context_gbr_peek_tab(struct isl_context *context)
2504 struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
2505 return cgbr->tab;
2508 /* Initialize the "shifted" tableau of the context, which
2509 * contains the constraints of the original tableau shifted
2510 * by the sum of all negative coefficients. This ensures
2511 * that any rational point in the shifted tableau can
2512 * be rounded up to yield an integer point in the original tableau.
2514 static void gbr_init_shifted(struct isl_context_gbr *cgbr)
2516 int i, j;
2517 struct isl_vec *cst;
2518 struct isl_basic_set *bset = isl_tab_peek_bset(cgbr->tab);
2519 unsigned dim = isl_basic_set_total_dim(bset);
2521 cst = isl_vec_alloc(cgbr->tab->mat->ctx, bset->n_ineq);
2522 if (!cst)
2523 return;
2525 for (i = 0; i < bset->n_ineq; ++i) {
2526 isl_int_set(cst->el[i], bset->ineq[i][0]);
2527 for (j = 0; j < dim; ++j) {
2528 if (!isl_int_is_neg(bset->ineq[i][1 + j]))
2529 continue;
2530 isl_int_add(bset->ineq[i][0], bset->ineq[i][0],
2531 bset->ineq[i][1 + j]);
2535 cgbr->shifted = isl_tab_from_basic_set(bset);
2537 for (i = 0; i < bset->n_ineq; ++i)
2538 isl_int_set(bset->ineq[i][0], cst->el[i]);
2540 isl_vec_free(cst);
2543 /* Check if the shifted tableau is non-empty, and if so
2544 * use the sample point to construct an integer point
2545 * of the context tableau.
2547 static struct isl_vec *gbr_get_shifted_sample(struct isl_context_gbr *cgbr)
2549 struct isl_vec *sample;
2551 if (!cgbr->shifted)
2552 gbr_init_shifted(cgbr);
2553 if (!cgbr->shifted)
2554 return NULL;
2555 if (cgbr->shifted->empty)
2556 return isl_vec_alloc(cgbr->tab->mat->ctx, 0);
2558 sample = isl_tab_get_sample_value(cgbr->shifted);
2559 sample = isl_vec_ceil(sample);
2561 return sample;
2564 static struct isl_basic_set *drop_constant_terms(struct isl_basic_set *bset)
2566 int i;
2568 if (!bset)
2569 return NULL;
2571 for (i = 0; i < bset->n_eq; ++i)
2572 isl_int_set_si(bset->eq[i][0], 0);
2574 for (i = 0; i < bset->n_ineq; ++i)
2575 isl_int_set_si(bset->ineq[i][0], 0);
2577 return bset;
2580 static int use_shifted(struct isl_context_gbr *cgbr)
2582 return cgbr->tab->bmap->n_eq == 0 && cgbr->tab->bmap->n_div == 0;
2585 static struct isl_vec *gbr_get_sample(struct isl_context_gbr *cgbr)
2587 struct isl_basic_set *bset;
2588 struct isl_basic_set *cone;
2590 if (isl_tab_sample_is_integer(cgbr->tab))
2591 return isl_tab_get_sample_value(cgbr->tab);
2593 if (use_shifted(cgbr)) {
2594 struct isl_vec *sample;
2596 sample = gbr_get_shifted_sample(cgbr);
2597 if (!sample || sample->size > 0)
2598 return sample;
2600 isl_vec_free(sample);
2603 if (!cgbr->cone) {
2604 bset = isl_tab_peek_bset(cgbr->tab);
2605 cgbr->cone = isl_tab_from_recession_cone(bset);
2606 if (!cgbr->cone)
2607 return NULL;
2608 if (isl_tab_track_bset(cgbr->cone, isl_basic_set_dup(bset)) < 0)
2609 return NULL;
2611 cgbr->cone = isl_tab_detect_implicit_equalities(cgbr->cone);
2612 if (!cgbr->cone)
2613 return NULL;
2615 if (cgbr->cone->n_dead == cgbr->cone->n_col) {
2616 struct isl_vec *sample;
2617 struct isl_tab_undo *snap;
2619 if (cgbr->tab->basis) {
2620 if (cgbr->tab->basis->n_col != 1 + cgbr->tab->n_var) {
2621 isl_mat_free(cgbr->tab->basis);
2622 cgbr->tab->basis = NULL;
2624 cgbr->tab->n_zero = 0;
2625 cgbr->tab->n_unbounded = 0;
2628 snap = isl_tab_snap(cgbr->tab);
2630 sample = isl_tab_sample(cgbr->tab);
2632 if (isl_tab_rollback(cgbr->tab, snap) < 0) {
2633 isl_vec_free(sample);
2634 return NULL;
2637 return sample;
2640 cone = isl_basic_set_dup(isl_tab_peek_bset(cgbr->cone));
2641 cone = drop_constant_terms(cone);
2642 cone = isl_basic_set_update_from_tab(cone, cgbr->cone);
2643 cone = isl_basic_set_underlying_set(cone);
2644 cone = isl_basic_set_gauss(cone, NULL);
2646 bset = isl_basic_set_dup(isl_tab_peek_bset(cgbr->tab));
2647 bset = isl_basic_set_update_from_tab(bset, cgbr->tab);
2648 bset = isl_basic_set_underlying_set(bset);
2649 bset = isl_basic_set_gauss(bset, NULL);
2651 return isl_basic_set_sample_with_cone(bset, cone);
2654 static void check_gbr_integer_feasible(struct isl_context_gbr *cgbr)
2656 struct isl_vec *sample;
2658 if (!cgbr->tab)
2659 return;
2661 if (cgbr->tab->empty)
2662 return;
2664 sample = gbr_get_sample(cgbr);
2665 if (!sample)
2666 goto error;
2668 if (sample->size == 0) {
2669 isl_vec_free(sample);
2670 if (isl_tab_mark_empty(cgbr->tab) < 0)
2671 goto error;
2672 return;
2675 cgbr->tab = isl_tab_add_sample(cgbr->tab, sample);
2677 return;
2678 error:
2679 isl_tab_free(cgbr->tab);
2680 cgbr->tab = NULL;
2683 static struct isl_tab *add_gbr_eq(struct isl_tab *tab, isl_int *eq)
2685 int r;
2687 if (!tab)
2688 return NULL;
2690 if (isl_tab_extend_cons(tab, 2) < 0)
2691 goto error;
2693 tab = isl_tab_add_eq(tab, eq);
2695 return tab;
2696 error:
2697 isl_tab_free(tab);
2698 return NULL;
2701 static void context_gbr_add_eq(struct isl_context *context, isl_int *eq,
2702 int check, int update)
2704 struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
2706 cgbr->tab = add_gbr_eq(cgbr->tab, eq);
2708 if (cgbr->cone && cgbr->cone->n_col != cgbr->cone->n_dead) {
2709 if (isl_tab_extend_cons(cgbr->cone, 2) < 0)
2710 goto error;
2711 cgbr->cone = isl_tab_add_eq(cgbr->cone, eq);
2714 if (check) {
2715 int v = tab_has_valid_sample(cgbr->tab, eq, 1);
2716 if (v < 0)
2717 goto error;
2718 if (!v)
2719 check_gbr_integer_feasible(cgbr);
2721 if (update)
2722 cgbr->tab = check_samples(cgbr->tab, eq, 1);
2723 return;
2724 error:
2725 isl_tab_free(cgbr->tab);
2726 cgbr->tab = NULL;
2729 static void add_gbr_ineq(struct isl_context_gbr *cgbr, isl_int *ineq)
2731 if (!cgbr->tab)
2732 return;
2734 if (isl_tab_extend_cons(cgbr->tab, 1) < 0)
2735 goto error;
2737 if (isl_tab_add_ineq(cgbr->tab, ineq) < 0)
2738 goto error;
2740 if (cgbr->shifted && !cgbr->shifted->empty && use_shifted(cgbr)) {
2741 int i;
2742 unsigned dim;
2743 dim = isl_basic_map_total_dim(cgbr->tab->bmap);
2745 if (isl_tab_extend_cons(cgbr->shifted, 1) < 0)
2746 goto error;
2748 for (i = 0; i < dim; ++i) {
2749 if (!isl_int_is_neg(ineq[1 + i]))
2750 continue;
2751 isl_int_add(ineq[0], ineq[0], ineq[1 + i]);
2754 if (isl_tab_add_ineq(cgbr->shifted, ineq) < 0)
2755 goto error;
2757 for (i = 0; i < dim; ++i) {
2758 if (!isl_int_is_neg(ineq[1 + i]))
2759 continue;
2760 isl_int_sub(ineq[0], ineq[0], ineq[1 + i]);
2764 if (cgbr->cone && cgbr->cone->n_col != cgbr->cone->n_dead) {
2765 if (isl_tab_extend_cons(cgbr->cone, 1) < 0)
2766 goto error;
2767 if (isl_tab_add_ineq(cgbr->cone, ineq) < 0)
2768 goto error;
2771 return;
2772 error:
2773 isl_tab_free(cgbr->tab);
2774 cgbr->tab = NULL;
2777 static void context_gbr_add_ineq(struct isl_context *context, isl_int *ineq,
2778 int check, int update)
2780 struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
2782 add_gbr_ineq(cgbr, ineq);
2783 if (!cgbr->tab)
2784 return;
2786 if (check) {
2787 int v = tab_has_valid_sample(cgbr->tab, ineq, 0);
2788 if (v < 0)
2789 goto error;
2790 if (!v)
2791 check_gbr_integer_feasible(cgbr);
2793 if (update)
2794 cgbr->tab = check_samples(cgbr->tab, ineq, 0);
2795 return;
2796 error:
2797 isl_tab_free(cgbr->tab);
2798 cgbr->tab = NULL;
2801 static int context_gbr_add_ineq_wrap(void *user, isl_int *ineq)
2803 struct isl_context *context = (struct isl_context *)user;
2804 context_gbr_add_ineq(context, ineq, 0, 0);
2805 return context->op->is_ok(context) ? 0 : -1;
2808 static enum isl_tab_row_sign context_gbr_ineq_sign(struct isl_context *context,
2809 isl_int *ineq, int strict)
2811 struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
2812 return tab_ineq_sign(cgbr->tab, ineq, strict);
2815 /* Check whether "ineq" can be added to the tableau without rendering
2816 * it infeasible.
2818 static int context_gbr_test_ineq(struct isl_context *context, isl_int *ineq)
2820 struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
2821 struct isl_tab_undo *snap;
2822 struct isl_tab_undo *shifted_snap = NULL;
2823 struct isl_tab_undo *cone_snap = NULL;
2824 int feasible;
2826 if (!cgbr->tab)
2827 return -1;
2829 if (isl_tab_extend_cons(cgbr->tab, 1) < 0)
2830 return -1;
2832 snap = isl_tab_snap(cgbr->tab);
2833 if (cgbr->shifted)
2834 shifted_snap = isl_tab_snap(cgbr->shifted);
2835 if (cgbr->cone)
2836 cone_snap = isl_tab_snap(cgbr->cone);
2837 add_gbr_ineq(cgbr, ineq);
2838 check_gbr_integer_feasible(cgbr);
2839 if (!cgbr->tab)
2840 return -1;
2841 feasible = !cgbr->tab->empty;
2842 if (isl_tab_rollback(cgbr->tab, snap) < 0)
2843 return -1;
2844 if (shifted_snap) {
2845 if (isl_tab_rollback(cgbr->shifted, shifted_snap))
2846 return -1;
2847 } else if (cgbr->shifted) {
2848 isl_tab_free(cgbr->shifted);
2849 cgbr->shifted = NULL;
2851 if (cone_snap) {
2852 if (isl_tab_rollback(cgbr->cone, cone_snap))
2853 return -1;
2854 } else if (cgbr->cone) {
2855 isl_tab_free(cgbr->cone);
2856 cgbr->cone = NULL;
2859 return feasible;
2862 /* Return the column of the last of the variables associated to
2863 * a column that has a non-zero coefficient.
2864 * This function is called in a context where only coefficients
2865 * of parameters or divs can be non-zero.
2867 static int last_non_zero_var_col(struct isl_tab *tab, isl_int *p)
2869 int i;
2870 int col;
2871 unsigned dim = tab->n_var - tab->n_param - tab->n_div;
2873 if (tab->n_var == 0)
2874 return -1;
2876 for (i = tab->n_var - 1; i >= 0; --i) {
2877 if (i >= tab->n_param && i < tab->n_var - tab->n_div)
2878 continue;
2879 if (tab->var[i].is_row)
2880 continue;
2881 col = tab->var[i].index;
2882 if (!isl_int_is_zero(p[col]))
2883 return col;
2886 return -1;
2889 /* Look through all the recently added equalities in the context
2890 * to see if we can propagate any of them to the main tableau.
2892 * The newly added equalities in the context are encoded as pairs
2893 * of inequalities starting at inequality "first".
2895 * We tentatively add each of these equalities to the main tableau
2896 * and if this happens to result in a row with a final coefficient
2897 * that is one or negative one, we use it to kill a column
2898 * in the main tableau. Otherwise, we discard the tentatively
2899 * added row.
2901 static void propagate_equalities(struct isl_context_gbr *cgbr,
2902 struct isl_tab *tab, unsigned first)
2904 int i;
2905 struct isl_vec *eq = NULL;
2907 eq = isl_vec_alloc(tab->mat->ctx, 1 + tab->n_var);
2908 if (!eq)
2909 goto error;
2911 if (isl_tab_extend_cons(tab, (cgbr->tab->bmap->n_ineq - first)/2) < 0)
2912 goto error;
2914 isl_seq_clr(eq->el + 1 + tab->n_param,
2915 tab->n_var - tab->n_param - tab->n_div);
2916 for (i = first; i < cgbr->tab->bmap->n_ineq; i += 2) {
2917 int j;
2918 int r;
2919 struct isl_tab_undo *snap;
2920 snap = isl_tab_snap(tab);
2922 isl_seq_cpy(eq->el, cgbr->tab->bmap->ineq[i], 1 + tab->n_param);
2923 isl_seq_cpy(eq->el + 1 + tab->n_var - tab->n_div,
2924 cgbr->tab->bmap->ineq[i] + 1 + tab->n_param,
2925 tab->n_div);
2927 r = isl_tab_add_row(tab, eq->el);
2928 if (r < 0)
2929 goto error;
2930 r = tab->con[r].index;
2931 j = last_non_zero_var_col(tab, tab->mat->row[r] + 2 + tab->M);
2932 if (j < 0 || j < tab->n_dead ||
2933 !isl_int_is_one(tab->mat->row[r][0]) ||
2934 (!isl_int_is_one(tab->mat->row[r][2 + tab->M + j]) &&
2935 !isl_int_is_negone(tab->mat->row[r][2 + tab->M + j]))) {
2936 if (isl_tab_rollback(tab, snap) < 0)
2937 goto error;
2938 continue;
2940 if (isl_tab_pivot(tab, r, j) < 0)
2941 goto error;
2942 if (isl_tab_kill_col(tab, j) < 0)
2943 goto error;
2945 tab = restore_lexmin(tab);
2948 isl_vec_free(eq);
2950 return;
2951 error:
2952 isl_vec_free(eq);
2953 isl_tab_free(cgbr->tab);
2954 cgbr->tab = NULL;
2957 static int context_gbr_detect_equalities(struct isl_context *context,
2958 struct isl_tab *tab)
2960 struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
2961 struct isl_ctx *ctx;
2962 int i;
2963 enum isl_lp_result res;
2964 unsigned n_ineq;
2966 ctx = cgbr->tab->mat->ctx;
2968 if (!cgbr->cone) {
2969 struct isl_basic_set *bset = isl_tab_peek_bset(cgbr->tab);
2970 cgbr->cone = isl_tab_from_recession_cone(bset);
2971 if (!cgbr->cone)
2972 goto error;
2973 if (isl_tab_track_bset(cgbr->cone, isl_basic_set_dup(bset)) < 0)
2974 goto error;
2976 cgbr->cone = isl_tab_detect_implicit_equalities(cgbr->cone);
2978 n_ineq = cgbr->tab->bmap->n_ineq;
2979 cgbr->tab = isl_tab_detect_equalities(cgbr->tab, cgbr->cone);
2980 if (cgbr->tab && cgbr->tab->bmap->n_ineq > n_ineq)
2981 propagate_equalities(cgbr, tab, n_ineq);
2983 return 0;
2984 error:
2985 isl_tab_free(cgbr->tab);
2986 cgbr->tab = NULL;
2987 return -1;
2990 static int context_gbr_get_div(struct isl_context *context, struct isl_tab *tab,
2991 struct isl_vec *div)
2993 return get_div(tab, context, div);
2996 static int context_gbr_add_div(struct isl_context *context, struct isl_vec *div)
2998 struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
2999 if (cgbr->cone) {
3000 int k;
3002 if (isl_tab_extend_cons(cgbr->cone, 3) < 0)
3003 return -1;
3004 if (isl_tab_extend_vars(cgbr->cone, 1) < 0)
3005 return -1;
3006 if (isl_tab_allocate_var(cgbr->cone) <0)
3007 return -1;
3009 cgbr->cone->bmap = isl_basic_map_extend_dim(cgbr->cone->bmap,
3010 isl_basic_map_get_dim(cgbr->cone->bmap), 1, 0, 2);
3011 k = isl_basic_map_alloc_div(cgbr->cone->bmap);
3012 if (k < 0)
3013 return -1;
3014 isl_seq_cpy(cgbr->cone->bmap->div[k], div->el, div->size);
3015 if (isl_tab_push(cgbr->cone, isl_tab_undo_bmap_div) < 0)
3016 return -1;
3018 return context_tab_add_div(cgbr->tab, div,
3019 context_gbr_add_ineq_wrap, context);
3022 static int context_gbr_best_split(struct isl_context *context,
3023 struct isl_tab *tab)
3025 struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3026 struct isl_tab_undo *snap;
3027 int r;
3029 snap = isl_tab_snap(cgbr->tab);
3030 r = best_split(tab, cgbr->tab);
3032 if (isl_tab_rollback(cgbr->tab, snap) < 0)
3033 return -1;
3035 return r;
3038 static int context_gbr_is_empty(struct isl_context *context)
3040 struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3041 if (!cgbr->tab)
3042 return -1;
3043 return cgbr->tab->empty;
3046 struct isl_gbr_tab_undo {
3047 struct isl_tab_undo *tab_snap;
3048 struct isl_tab_undo *shifted_snap;
3049 struct isl_tab_undo *cone_snap;
3052 static void *context_gbr_save(struct isl_context *context)
3054 struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3055 struct isl_gbr_tab_undo *snap;
3057 snap = isl_alloc_type(cgbr->tab->mat->ctx, struct isl_gbr_tab_undo);
3058 if (!snap)
3059 return NULL;
3061 snap->tab_snap = isl_tab_snap(cgbr->tab);
3062 if (isl_tab_save_samples(cgbr->tab) < 0)
3063 goto error;
3065 if (cgbr->shifted)
3066 snap->shifted_snap = isl_tab_snap(cgbr->shifted);
3067 else
3068 snap->shifted_snap = NULL;
3070 if (cgbr->cone)
3071 snap->cone_snap = isl_tab_snap(cgbr->cone);
3072 else
3073 snap->cone_snap = NULL;
3075 return snap;
3076 error:
3077 free(snap);
3078 return NULL;
3081 static void context_gbr_restore(struct isl_context *context, void *save)
3083 struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3084 struct isl_gbr_tab_undo *snap = (struct isl_gbr_tab_undo *)save;
3085 if (!snap)
3086 goto error;
3087 if (isl_tab_rollback(cgbr->tab, snap->tab_snap) < 0) {
3088 isl_tab_free(cgbr->tab);
3089 cgbr->tab = NULL;
3092 if (snap->shifted_snap) {
3093 if (isl_tab_rollback(cgbr->shifted, snap->shifted_snap) < 0)
3094 goto error;
3095 } else if (cgbr->shifted) {
3096 isl_tab_free(cgbr->shifted);
3097 cgbr->shifted = NULL;
3100 if (snap->cone_snap) {
3101 if (isl_tab_rollback(cgbr->cone, snap->cone_snap) < 0)
3102 goto error;
3103 } else if (cgbr->cone) {
3104 isl_tab_free(cgbr->cone);
3105 cgbr->cone = NULL;
3108 free(snap);
3110 return;
3111 error:
3112 free(snap);
3113 isl_tab_free(cgbr->tab);
3114 cgbr->tab = NULL;
3117 static int context_gbr_is_ok(struct isl_context *context)
3119 struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3120 return !!cgbr->tab;
3123 static void context_gbr_invalidate(struct isl_context *context)
3125 struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3126 isl_tab_free(cgbr->tab);
3127 cgbr->tab = NULL;
3130 static void context_gbr_free(struct isl_context *context)
3132 struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3133 isl_tab_free(cgbr->tab);
3134 isl_tab_free(cgbr->shifted);
3135 isl_tab_free(cgbr->cone);
3136 free(cgbr);
3139 struct isl_context_op isl_context_gbr_op = {
3140 context_gbr_detect_nonnegative_parameters,
3141 context_gbr_peek_basic_set,
3142 context_gbr_peek_tab,
3143 context_gbr_add_eq,
3144 context_gbr_add_ineq,
3145 context_gbr_ineq_sign,
3146 context_gbr_test_ineq,
3147 context_gbr_get_div,
3148 context_gbr_add_div,
3149 context_gbr_detect_equalities,
3150 context_gbr_best_split,
3151 context_gbr_is_empty,
3152 context_gbr_is_ok,
3153 context_gbr_save,
3154 context_gbr_restore,
3155 context_gbr_invalidate,
3156 context_gbr_free,
3159 static struct isl_context *isl_context_gbr_alloc(struct isl_basic_set *dom)
3161 struct isl_context_gbr *cgbr;
3163 if (!dom)
3164 return NULL;
3166 cgbr = isl_calloc_type(dom->ctx, struct isl_context_gbr);
3167 if (!cgbr)
3168 return NULL;
3170 cgbr->context.op = &isl_context_gbr_op;
3172 cgbr->shifted = NULL;
3173 cgbr->cone = NULL;
3174 cgbr->tab = isl_tab_from_basic_set(dom);
3175 cgbr->tab = isl_tab_init_samples(cgbr->tab);
3176 if (!cgbr->tab)
3177 goto error;
3178 if (isl_tab_track_bset(cgbr->tab,
3179 isl_basic_set_cow(isl_basic_set_copy(dom))) < 0)
3180 goto error;
3181 check_gbr_integer_feasible(cgbr);
3183 return &cgbr->context;
3184 error:
3185 cgbr->context.op->free(&cgbr->context);
3186 return NULL;
3189 static struct isl_context *isl_context_alloc(struct isl_basic_set *dom)
3191 if (!dom)
3192 return NULL;
3194 if (dom->ctx->opt->context == ISL_CONTEXT_LEXMIN)
3195 return isl_context_lex_alloc(dom);
3196 else
3197 return isl_context_gbr_alloc(dom);
3200 /* Construct an isl_sol_map structure for accumulating the solution.
3201 * If track_empty is set, then we also keep track of the parts
3202 * of the context where there is no solution.
3203 * If max is set, then we are solving a maximization, rather than
3204 * a minimization problem, which means that the variables in the
3205 * tableau have value "M - x" rather than "M + x".
3207 static struct isl_sol_map *sol_map_init(struct isl_basic_map *bmap,
3208 struct isl_basic_set *dom, int track_empty, int max)
3210 struct isl_sol_map *sol_map;
3212 sol_map = isl_calloc_type(bset->ctx, struct isl_sol_map);
3213 if (!sol_map)
3214 goto error;
3216 sol_map->sol.rational = ISL_F_ISSET(bmap, ISL_BASIC_MAP_RATIONAL);
3217 sol_map->sol.dec_level.callback.run = &sol_dec_level_wrap;
3218 sol_map->sol.dec_level.sol = &sol_map->sol;
3219 sol_map->sol.max = max;
3220 sol_map->sol.n_out = isl_basic_map_dim(bmap, isl_dim_out);
3221 sol_map->sol.add = &sol_map_add_wrap;
3222 sol_map->sol.add_empty = track_empty ? &sol_map_add_empty_wrap : NULL;
3223 sol_map->sol.free = &sol_map_free_wrap;
3224 sol_map->map = isl_map_alloc_dim(isl_basic_map_get_dim(bmap), 1,
3225 ISL_MAP_DISJOINT);
3226 if (!sol_map->map)
3227 goto error;
3229 sol_map->sol.context = isl_context_alloc(dom);
3230 if (!sol_map->sol.context)
3231 goto error;
3233 if (track_empty) {
3234 sol_map->empty = isl_set_alloc_dim(isl_basic_set_get_dim(dom),
3235 1, ISL_SET_DISJOINT);
3236 if (!sol_map->empty)
3237 goto error;
3240 isl_basic_set_free(dom);
3241 return sol_map;
3242 error:
3243 isl_basic_set_free(dom);
3244 sol_map_free(sol_map);
3245 return NULL;
3248 /* Check whether all coefficients of (non-parameter) variables
3249 * are non-positive, meaning that no pivots can be performed on the row.
3251 static int is_critical(struct isl_tab *tab, int row)
3253 int j;
3254 unsigned off = 2 + tab->M;
3256 for (j = tab->n_dead; j < tab->n_col; ++j) {
3257 if (tab->col_var[j] >= 0 &&
3258 (tab->col_var[j] < tab->n_param ||
3259 tab->col_var[j] >= tab->n_var - tab->n_div))
3260 continue;
3262 if (isl_int_is_pos(tab->mat->row[row][off + j]))
3263 return 0;
3266 return 1;
3269 /* Check whether the inequality represented by vec is strict over the integers,
3270 * i.e., there are no integer values satisfying the constraint with
3271 * equality. This happens if the gcd of the coefficients is not a divisor
3272 * of the constant term. If so, scale the constraint down by the gcd
3273 * of the coefficients.
3275 static int is_strict(struct isl_vec *vec)
3277 isl_int gcd;
3278 int strict = 0;
3280 isl_int_init(gcd);
3281 isl_seq_gcd(vec->el + 1, vec->size - 1, &gcd);
3282 if (!isl_int_is_one(gcd)) {
3283 strict = !isl_int_is_divisible_by(vec->el[0], gcd);
3284 isl_int_fdiv_q(vec->el[0], vec->el[0], gcd);
3285 isl_seq_scale_down(vec->el + 1, vec->el + 1, gcd, vec->size-1);
3287 isl_int_clear(gcd);
3289 return strict;
3292 /* Determine the sign of the given row of the main tableau.
3293 * The result is one of
3294 * isl_tab_row_pos: always non-negative; no pivot needed
3295 * isl_tab_row_neg: always non-positive; pivot
3296 * isl_tab_row_any: can be both positive and negative; split
3298 * We first handle some simple cases
3299 * - the row sign may be known already
3300 * - the row may be obviously non-negative
3301 * - the parametric constant may be equal to that of another row
3302 * for which we know the sign. This sign will be either "pos" or
3303 * "any". If it had been "neg" then we would have pivoted before.
3305 * If none of these cases hold, we check the value of the row for each
3306 * of the currently active samples. Based on the signs of these values
3307 * we make an initial determination of the sign of the row.
3309 * all zero -> unk(nown)
3310 * all non-negative -> pos
3311 * all non-positive -> neg
3312 * both negative and positive -> all
3314 * If we end up with "all", we are done.
3315 * Otherwise, we perform a check for positive and/or negative
3316 * values as follows.
3318 * samples neg unk pos
3319 * <0 ? Y N Y N
3320 * pos any pos
3321 * >0 ? Y N Y N
3322 * any neg any neg
3324 * There is no special sign for "zero", because we can usually treat zero
3325 * as either non-negative or non-positive, whatever works out best.
3326 * However, if the row is "critical", meaning that pivoting is impossible
3327 * then we don't want to limp zero with the non-positive case, because
3328 * then we we would lose the solution for those values of the parameters
3329 * where the value of the row is zero. Instead, we treat 0 as non-negative
3330 * ensuring a split if the row can attain both zero and negative values.
3331 * The same happens when the original constraint was one that could not
3332 * be satisfied with equality by any integer values of the parameters.
3333 * In this case, we normalize the constraint, but then a value of zero
3334 * for the normalized constraint is actually a positive value for the
3335 * original constraint, so again we need to treat zero as non-negative.
3336 * In both these cases, we have the following decision tree instead:
3338 * all non-negative -> pos
3339 * all negative -> neg
3340 * both negative and non-negative -> all
3342 * samples neg pos
3343 * <0 ? Y N
3344 * any pos
3345 * >=0 ? Y N
3346 * any neg
3348 static enum isl_tab_row_sign row_sign(struct isl_tab *tab,
3349 struct isl_sol *sol, int row)
3351 struct isl_vec *ineq = NULL;
3352 int res = isl_tab_row_unknown;
3353 int critical;
3354 int strict;
3355 int row2;
3357 if (tab->row_sign[row] != isl_tab_row_unknown)
3358 return tab->row_sign[row];
3359 if (is_obviously_nonneg(tab, row))
3360 return isl_tab_row_pos;
3361 for (row2 = tab->n_redundant; row2 < tab->n_row; ++row2) {
3362 if (tab->row_sign[row2] == isl_tab_row_unknown)
3363 continue;
3364 if (identical_parameter_line(tab, row, row2))
3365 return tab->row_sign[row2];
3368 critical = is_critical(tab, row);
3370 ineq = get_row_parameter_ineq(tab, row);
3371 if (!ineq)
3372 goto error;
3374 strict = is_strict(ineq);
3376 res = sol->context->op->ineq_sign(sol->context, ineq->el,
3377 critical || strict);
3379 if (res == isl_tab_row_unknown || res == isl_tab_row_pos) {
3380 /* test for negative values */
3381 int feasible;
3382 isl_seq_neg(ineq->el, ineq->el, ineq->size);
3383 isl_int_sub_ui(ineq->el[0], ineq->el[0], 1);
3385 feasible = sol->context->op->test_ineq(sol->context, ineq->el);
3386 if (feasible < 0)
3387 goto error;
3388 if (!feasible)
3389 res = isl_tab_row_pos;
3390 else
3391 res = (res == isl_tab_row_unknown) ? isl_tab_row_neg
3392 : isl_tab_row_any;
3393 if (res == isl_tab_row_neg) {
3394 isl_seq_neg(ineq->el, ineq->el, ineq->size);
3395 isl_int_sub_ui(ineq->el[0], ineq->el[0], 1);
3399 if (res == isl_tab_row_neg) {
3400 /* test for positive values */
3401 int feasible;
3402 if (!critical && !strict)
3403 isl_int_sub_ui(ineq->el[0], ineq->el[0], 1);
3405 feasible = sol->context->op->test_ineq(sol->context, ineq->el);
3406 if (feasible < 0)
3407 goto error;
3408 if (feasible)
3409 res = isl_tab_row_any;
3412 isl_vec_free(ineq);
3413 return res;
3414 error:
3415 isl_vec_free(ineq);
3416 return 0;
3419 static void find_solutions(struct isl_sol *sol, struct isl_tab *tab);
3421 /* Find solutions for values of the parameters that satisfy the given
3422 * inequality.
3424 * We currently take a snapshot of the context tableau that is reset
3425 * when we return from this function, while we make a copy of the main
3426 * tableau, leaving the original main tableau untouched.
3427 * These are fairly arbitrary choices. Making a copy also of the context
3428 * tableau would obviate the need to undo any changes made to it later,
3429 * while taking a snapshot of the main tableau could reduce memory usage.
3430 * If we were to switch to taking a snapshot of the main tableau,
3431 * we would have to keep in mind that we need to save the row signs
3432 * and that we need to do this before saving the current basis
3433 * such that the basis has been restore before we restore the row signs.
3435 static void find_in_pos(struct isl_sol *sol, struct isl_tab *tab, isl_int *ineq)
3437 void *saved;
3439 if (!sol->context)
3440 goto error;
3441 saved = sol->context->op->save(sol->context);
3443 tab = isl_tab_dup(tab);
3444 if (!tab)
3445 goto error;
3447 sol->context->op->add_ineq(sol->context, ineq, 0, 1);
3449 find_solutions(sol, tab);
3451 sol->context->op->restore(sol->context, saved);
3452 return;
3453 error:
3454 sol->error = 1;
3457 /* Record the absence of solutions for those values of the parameters
3458 * that do not satisfy the given inequality with equality.
3460 static void no_sol_in_strict(struct isl_sol *sol,
3461 struct isl_tab *tab, struct isl_vec *ineq)
3463 int empty;
3464 void *saved;
3466 if (!sol->context)
3467 goto error;
3468 saved = sol->context->op->save(sol->context);
3470 isl_int_sub_ui(ineq->el[0], ineq->el[0], 1);
3472 sol->context->op->add_ineq(sol->context, ineq->el, 1, 0);
3473 if (!sol->context)
3474 goto error;
3476 empty = tab->empty;
3477 tab->empty = 1;
3478 sol_add(sol, tab);
3479 tab->empty = empty;
3481 isl_int_add_ui(ineq->el[0], ineq->el[0], 1);
3483 sol->context->op->restore(sol->context, saved);
3484 return;
3485 error:
3486 sol->error = 1;
3489 /* Compute the lexicographic minimum of the set represented by the main
3490 * tableau "tab" within the context "sol->context_tab".
3491 * On entry the sample value of the main tableau is lexicographically
3492 * less than or equal to this lexicographic minimum.
3493 * Pivots are performed until a feasible point is found, which is then
3494 * necessarily equal to the minimum, or until the tableau is found to
3495 * be infeasible. Some pivots may need to be performed for only some
3496 * feasible values of the context tableau. If so, the context tableau
3497 * is split into a part where the pivot is needed and a part where it is not.
3499 * Whenever we enter the main loop, the main tableau is such that no
3500 * "obvious" pivots need to be performed on it, where "obvious" means
3501 * that the given row can be seen to be negative without looking at
3502 * the context tableau. In particular, for non-parametric problems,
3503 * no pivots need to be performed on the main tableau.
3504 * The caller of find_solutions is responsible for making this property
3505 * hold prior to the first iteration of the loop, while restore_lexmin
3506 * is called before every other iteration.
3508 * Inside the main loop, we first examine the signs of the rows of
3509 * the main tableau within the context of the context tableau.
3510 * If we find a row that is always non-positive for all values of
3511 * the parameters satisfying the context tableau and negative for at
3512 * least one value of the parameters, we perform the appropriate pivot
3513 * and start over. An exception is the case where no pivot can be
3514 * performed on the row. In this case, we require that the sign of
3515 * the row is negative for all values of the parameters (rather than just
3516 * non-positive). This special case is handled inside row_sign, which
3517 * will say that the row can have any sign if it determines that it can
3518 * attain both negative and zero values.
3520 * If we can't find a row that always requires a pivot, but we can find
3521 * one or more rows that require a pivot for some values of the parameters
3522 * (i.e., the row can attain both positive and negative signs), then we split
3523 * the context tableau into two parts, one where we force the sign to be
3524 * non-negative and one where we force is to be negative.
3525 * The non-negative part is handled by a recursive call (through find_in_pos).
3526 * Upon returning from this call, we continue with the negative part and
3527 * perform the required pivot.
3529 * If no such rows can be found, all rows are non-negative and we have
3530 * found a (rational) feasible point. If we only wanted a rational point
3531 * then we are done.
3532 * Otherwise, we check if all values of the sample point of the tableau
3533 * are integral for the variables. If so, we have found the minimal
3534 * integral point and we are done.
3535 * If the sample point is not integral, then we need to make a distinction
3536 * based on whether the constant term is non-integral or the coefficients
3537 * of the parameters. Furthermore, in order to decide how to handle
3538 * the non-integrality, we also need to know whether the coefficients
3539 * of the other columns in the tableau are integral. This leads
3540 * to the following table. The first two rows do not correspond
3541 * to a non-integral sample point and are only mentioned for completeness.
3543 * constant parameters other
3545 * int int int |
3546 * int int rat | -> no problem
3548 * rat int int -> fail
3550 * rat int rat -> cut
3552 * int rat rat |
3553 * rat rat rat | -> parametric cut
3555 * int rat int |
3556 * rat rat int | -> split context
3558 * If the parametric constant is completely integral, then there is nothing
3559 * to be done. If the constant term is non-integral, but all the other
3560 * coefficient are integral, then there is nothing that can be done
3561 * and the tableau has no integral solution.
3562 * If, on the other hand, one or more of the other columns have rational
3563 * coeffcients, but the parameter coefficients are all integral, then
3564 * we can perform a regular (non-parametric) cut.
3565 * Finally, if there is any parameter coefficient that is non-integral,
3566 * then we need to involve the context tableau. There are two cases here.
3567 * If at least one other column has a rational coefficient, then we
3568 * can perform a parametric cut in the main tableau by adding a new
3569 * integer division in the context tableau.
3570 * If all other columns have integral coefficients, then we need to
3571 * enforce that the rational combination of parameters (c + \sum a_i y_i)/m
3572 * is always integral. We do this by introducing an integer division
3573 * q = floor((c + \sum a_i y_i)/m) and stipulating that its argument should
3574 * always be integral in the context tableau, i.e., m q = c + \sum a_i y_i.
3575 * Since q is expressed in the tableau as
3576 * c + \sum a_i y_i - m q >= 0
3577 * -c - \sum a_i y_i + m q + m - 1 >= 0
3578 * it is sufficient to add the inequality
3579 * -c - \sum a_i y_i + m q >= 0
3580 * In the part of the context where this inequality does not hold, the
3581 * main tableau is marked as being empty.
3583 static void find_solutions(struct isl_sol *sol, struct isl_tab *tab)
3585 struct isl_context *context;
3587 if (!tab || sol->error)
3588 goto error;
3590 context = sol->context;
3592 if (tab->empty)
3593 goto done;
3594 if (context->op->is_empty(context))
3595 goto done;
3597 for (; tab && !tab->empty; tab = restore_lexmin(tab)) {
3598 int flags;
3599 int row;
3600 int sgn;
3601 int split = -1;
3602 int n_split = 0;
3604 for (row = tab->n_redundant; row < tab->n_row; ++row) {
3605 if (!isl_tab_var_from_row(tab, row)->is_nonneg)
3606 continue;
3607 sgn = row_sign(tab, sol, row);
3608 if (!sgn)
3609 goto error;
3610 tab->row_sign[row] = sgn;
3611 if (sgn == isl_tab_row_any)
3612 n_split++;
3613 if (sgn == isl_tab_row_any && split == -1)
3614 split = row;
3615 if (sgn == isl_tab_row_neg)
3616 break;
3618 if (row < tab->n_row)
3619 continue;
3620 if (split != -1) {
3621 struct isl_vec *ineq;
3622 if (n_split != 1)
3623 split = context->op->best_split(context, tab);
3624 if (split < 0)
3625 goto error;
3626 ineq = get_row_parameter_ineq(tab, split);
3627 if (!ineq)
3628 goto error;
3629 is_strict(ineq);
3630 for (row = tab->n_redundant; row < tab->n_row; ++row) {
3631 if (!isl_tab_var_from_row(tab, row)->is_nonneg)
3632 continue;
3633 if (tab->row_sign[row] == isl_tab_row_any)
3634 tab->row_sign[row] = isl_tab_row_unknown;
3636 tab->row_sign[split] = isl_tab_row_pos;
3637 sol_inc_level(sol);
3638 find_in_pos(sol, tab, ineq->el);
3639 tab->row_sign[split] = isl_tab_row_neg;
3640 row = split;
3641 isl_seq_neg(ineq->el, ineq->el, ineq->size);
3642 isl_int_sub_ui(ineq->el[0], ineq->el[0], 1);
3643 context->op->add_ineq(context, ineq->el, 0, 1);
3644 isl_vec_free(ineq);
3645 if (sol->error)
3646 goto error;
3647 continue;
3649 if (tab->rational)
3650 break;
3651 row = first_non_integer_row(tab, &flags);
3652 if (row < 0)
3653 break;
3654 if (ISL_FL_ISSET(flags, I_PAR)) {
3655 if (ISL_FL_ISSET(flags, I_VAR)) {
3656 if (isl_tab_mark_empty(tab) < 0)
3657 goto error;
3658 break;
3660 row = add_cut(tab, row);
3661 } else if (ISL_FL_ISSET(flags, I_VAR)) {
3662 struct isl_vec *div;
3663 struct isl_vec *ineq;
3664 int d;
3665 div = get_row_split_div(tab, row);
3666 if (!div)
3667 goto error;
3668 d = context->op->get_div(context, tab, div);
3669 isl_vec_free(div);
3670 if (d < 0)
3671 goto error;
3672 ineq = ineq_for_div(context->op->peek_basic_set(context), d);
3673 sol_inc_level(sol);
3674 no_sol_in_strict(sol, tab, ineq);
3675 isl_seq_neg(ineq->el, ineq->el, ineq->size);
3676 context->op->add_ineq(context, ineq->el, 1, 1);
3677 isl_vec_free(ineq);
3678 if (sol->error || !context->op->is_ok(context))
3679 goto error;
3680 tab = set_row_cst_to_div(tab, row, d);
3681 if (context->op->is_empty(context))
3682 break;
3683 } else
3684 row = add_parametric_cut(tab, row, context);
3685 if (row < 0)
3686 goto error;
3688 done:
3689 sol_add(sol, tab);
3690 isl_tab_free(tab);
3691 return;
3692 error:
3693 isl_tab_free(tab);
3694 sol_free(sol);
3697 /* Compute the lexicographic minimum of the set represented by the main
3698 * tableau "tab" within the context "sol->context_tab".
3700 * As a preprocessing step, we first transfer all the purely parametric
3701 * equalities from the main tableau to the context tableau, i.e.,
3702 * parameters that have been pivoted to a row.
3703 * These equalities are ignored by the main algorithm, because the
3704 * corresponding rows may not be marked as being non-negative.
3705 * In parts of the context where the added equality does not hold,
3706 * the main tableau is marked as being empty.
3708 static void find_solutions_main(struct isl_sol *sol, struct isl_tab *tab)
3710 int row;
3712 sol->level = 0;
3714 for (row = tab->n_redundant; row < tab->n_row; ++row) {
3715 int p;
3716 struct isl_vec *eq;
3718 if (tab->row_var[row] < 0)
3719 continue;
3720 if (tab->row_var[row] >= tab->n_param &&
3721 tab->row_var[row] < tab->n_var - tab->n_div)
3722 continue;
3723 if (tab->row_var[row] < tab->n_param)
3724 p = tab->row_var[row];
3725 else
3726 p = tab->row_var[row]
3727 + tab->n_param - (tab->n_var - tab->n_div);
3729 eq = isl_vec_alloc(tab->mat->ctx, 1+tab->n_param+tab->n_div);
3730 get_row_parameter_line(tab, row, eq->el);
3731 isl_int_neg(eq->el[1 + p], tab->mat->row[row][0]);
3732 eq = isl_vec_normalize(eq);
3734 sol_inc_level(sol);
3735 no_sol_in_strict(sol, tab, eq);
3737 isl_seq_neg(eq->el, eq->el, eq->size);
3738 sol_inc_level(sol);
3739 no_sol_in_strict(sol, tab, eq);
3740 isl_seq_neg(eq->el, eq->el, eq->size);
3742 sol->context->op->add_eq(sol->context, eq->el, 1, 1);
3744 isl_vec_free(eq);
3746 if (isl_tab_mark_redundant(tab, row) < 0)
3747 goto error;
3749 if (sol->context->op->is_empty(sol->context))
3750 break;
3752 row = tab->n_redundant - 1;
3755 find_solutions(sol, tab);
3757 sol->level = 0;
3758 sol_pop(sol);
3760 return;
3761 error:
3762 isl_tab_free(tab);
3763 sol_free(sol);
3766 static void sol_map_find_solutions(struct isl_sol_map *sol_map,
3767 struct isl_tab *tab)
3769 find_solutions_main(&sol_map->sol, tab);
3772 /* Check if integer division "div" of "dom" also occurs in "bmap".
3773 * If so, return its position within the divs.
3774 * If not, return -1.
3776 static int find_context_div(struct isl_basic_map *bmap,
3777 struct isl_basic_set *dom, unsigned div)
3779 int i;
3780 unsigned b_dim = isl_dim_total(bmap->dim);
3781 unsigned d_dim = isl_dim_total(dom->dim);
3783 if (isl_int_is_zero(dom->div[div][0]))
3784 return -1;
3785 if (isl_seq_first_non_zero(dom->div[div] + 2 + d_dim, dom->n_div) != -1)
3786 return -1;
3788 for (i = 0; i < bmap->n_div; ++i) {
3789 if (isl_int_is_zero(bmap->div[i][0]))
3790 continue;
3791 if (isl_seq_first_non_zero(bmap->div[i] + 2 + d_dim,
3792 (b_dim - d_dim) + bmap->n_div) != -1)
3793 continue;
3794 if (isl_seq_eq(bmap->div[i], dom->div[div], 2 + d_dim))
3795 return i;
3797 return -1;
3800 /* The correspondence between the variables in the main tableau,
3801 * the context tableau, and the input map and domain is as follows.
3802 * The first n_param and the last n_div variables of the main tableau
3803 * form the variables of the context tableau.
3804 * In the basic map, these n_param variables correspond to the
3805 * parameters and the input dimensions. In the domain, they correspond
3806 * to the parameters and the set dimensions.
3807 * The n_div variables correspond to the integer divisions in the domain.
3808 * To ensure that everything lines up, we may need to copy some of the
3809 * integer divisions of the domain to the map. These have to be placed
3810 * in the same order as those in the context and they have to be placed
3811 * after any other integer divisions that the map may have.
3812 * This function performs the required reordering.
3814 static struct isl_basic_map *align_context_divs(struct isl_basic_map *bmap,
3815 struct isl_basic_set *dom)
3817 int i;
3818 int common = 0;
3819 int other;
3821 for (i = 0; i < dom->n_div; ++i)
3822 if (find_context_div(bmap, dom, i) != -1)
3823 common++;
3824 other = bmap->n_div - common;
3825 if (dom->n_div - common > 0) {
3826 bmap = isl_basic_map_extend_dim(bmap, isl_dim_copy(bmap->dim),
3827 dom->n_div - common, 0, 0);
3828 if (!bmap)
3829 return NULL;
3831 for (i = 0; i < dom->n_div; ++i) {
3832 int pos = find_context_div(bmap, dom, i);
3833 if (pos < 0) {
3834 pos = isl_basic_map_alloc_div(bmap);
3835 if (pos < 0)
3836 goto error;
3837 isl_int_set_si(bmap->div[pos][0], 0);
3839 if (pos != other + i)
3840 isl_basic_map_swap_div(bmap, pos, other + i);
3842 return bmap;
3843 error:
3844 isl_basic_map_free(bmap);
3845 return NULL;
3848 /* Compute the lexicographic minimum (or maximum if "max" is set)
3849 * of "bmap" over the domain "dom" and return the result as a map.
3850 * If "empty" is not NULL, then *empty is assigned a set that
3851 * contains those parts of the domain where there is no solution.
3852 * If "bmap" is marked as rational (ISL_BASIC_MAP_RATIONAL),
3853 * then we compute the rational optimum. Otherwise, we compute
3854 * the integral optimum.
3856 * We perform some preprocessing. As the PILP solver does not
3857 * handle implicit equalities very well, we first make sure all
3858 * the equalities are explicitly available.
3859 * We also make sure the divs in the domain are properly order,
3860 * because they will be added one by one in the given order
3861 * during the construction of the solution map.
3863 struct isl_map *isl_tab_basic_map_partial_lexopt(
3864 struct isl_basic_map *bmap, struct isl_basic_set *dom,
3865 struct isl_set **empty, int max)
3867 struct isl_tab *tab;
3868 struct isl_map *result = NULL;
3869 struct isl_sol_map *sol_map = NULL;
3870 struct isl_context *context;
3871 struct isl_basic_map *eq;
3873 if (empty)
3874 *empty = NULL;
3875 if (!bmap || !dom)
3876 goto error;
3878 isl_assert(bmap->ctx,
3879 isl_basic_map_compatible_domain(bmap, dom), goto error);
3881 eq = isl_basic_map_copy(bmap);
3882 eq = isl_basic_map_intersect_domain(eq, isl_basic_set_copy(dom));
3883 eq = isl_basic_map_affine_hull(eq);
3884 bmap = isl_basic_map_intersect(bmap, eq);
3886 if (dom->n_div) {
3887 dom = isl_basic_set_order_divs(dom);
3888 bmap = align_context_divs(bmap, dom);
3890 sol_map = sol_map_init(bmap, dom, !!empty, max);
3891 if (!sol_map)
3892 goto error;
3894 context = sol_map->sol.context;
3895 if (isl_basic_set_fast_is_empty(context->op->peek_basic_set(context)))
3896 /* nothing */;
3897 else if (isl_basic_map_fast_is_empty(bmap))
3898 sol_map_add_empty_if_needed(sol_map,
3899 isl_basic_set_copy(context->op->peek_basic_set(context)));
3900 else {
3901 tab = tab_for_lexmin(bmap,
3902 context->op->peek_basic_set(context), 1, max);
3903 tab = context->op->detect_nonnegative_parameters(context, tab);
3904 sol_map_find_solutions(sol_map, tab);
3906 if (sol_map->sol.error)
3907 goto error;
3909 result = isl_map_copy(sol_map->map);
3910 if (empty)
3911 *empty = isl_set_copy(sol_map->empty);
3912 sol_free(&sol_map->sol);
3913 isl_basic_map_free(bmap);
3914 return result;
3915 error:
3916 sol_free(&sol_map->sol);
3917 isl_basic_map_free(bmap);
3918 return NULL;
3921 struct isl_sol_for {
3922 struct isl_sol sol;
3923 int (*fn)(__isl_take isl_basic_set *dom,
3924 __isl_take isl_mat *map, void *user);
3925 void *user;
3928 static void sol_for_free(struct isl_sol_for *sol_for)
3930 if (sol_for->sol.context)
3931 sol_for->sol.context->op->free(sol_for->sol.context);
3932 free(sol_for);
3935 static void sol_for_free_wrap(struct isl_sol *sol)
3937 sol_for_free((struct isl_sol_for *)sol);
3940 /* Add the solution identified by the tableau and the context tableau.
3942 * See documentation of sol_add for more details.
3944 * Instead of constructing a basic map, this function calls a user
3945 * defined function with the current context as a basic set and
3946 * an affine matrix reprenting the relation between the input and output.
3947 * The number of rows in this matrix is equal to one plus the number
3948 * of output variables. The number of columns is equal to one plus
3949 * the total dimension of the context, i.e., the number of parameters,
3950 * input variables and divs. Since some of the columns in the matrix
3951 * may refer to the divs, the basic set is not simplified.
3952 * (Simplification may reorder or remove divs.)
3954 static void sol_for_add(struct isl_sol_for *sol,
3955 struct isl_basic_set *dom, struct isl_mat *M)
3957 if (sol->sol.error || !dom || !M)
3958 goto error;
3960 dom = isl_basic_set_simplify(dom);
3961 dom = isl_basic_set_finalize(dom);
3963 if (sol->fn(isl_basic_set_copy(dom), isl_mat_copy(M), sol->user) < 0)
3964 goto error;
3966 isl_basic_set_free(dom);
3967 isl_mat_free(M);
3968 return;
3969 error:
3970 isl_basic_set_free(dom);
3971 isl_mat_free(M);
3972 sol->sol.error = 1;
3975 static void sol_for_add_wrap(struct isl_sol *sol,
3976 struct isl_basic_set *dom, struct isl_mat *M)
3978 sol_for_add((struct isl_sol_for *)sol, dom, M);
3981 static struct isl_sol_for *sol_for_init(struct isl_basic_map *bmap, int max,
3982 int (*fn)(__isl_take isl_basic_set *dom, __isl_take isl_mat *map,
3983 void *user),
3984 void *user)
3986 struct isl_sol_for *sol_for = NULL;
3987 struct isl_dim *dom_dim;
3988 struct isl_basic_set *dom = NULL;
3990 sol_for = isl_calloc_type(bset->ctx, struct isl_sol_for);
3991 if (!sol_for)
3992 goto error;
3994 dom_dim = isl_dim_domain(isl_dim_copy(bmap->dim));
3995 dom = isl_basic_set_universe(dom_dim);
3997 sol_for->sol.rational = ISL_F_ISSET(bmap, ISL_BASIC_MAP_RATIONAL);
3998 sol_for->sol.dec_level.callback.run = &sol_dec_level_wrap;
3999 sol_for->sol.dec_level.sol = &sol_for->sol;
4000 sol_for->fn = fn;
4001 sol_for->user = user;
4002 sol_for->sol.max = max;
4003 sol_for->sol.n_out = isl_basic_map_dim(bmap, isl_dim_out);
4004 sol_for->sol.add = &sol_for_add_wrap;
4005 sol_for->sol.add_empty = NULL;
4006 sol_for->sol.free = &sol_for_free_wrap;
4008 sol_for->sol.context = isl_context_alloc(dom);
4009 if (!sol_for->sol.context)
4010 goto error;
4012 isl_basic_set_free(dom);
4013 return sol_for;
4014 error:
4015 isl_basic_set_free(dom);
4016 sol_for_free(sol_for);
4017 return NULL;
4020 static void sol_for_find_solutions(struct isl_sol_for *sol_for,
4021 struct isl_tab *tab)
4023 find_solutions_main(&sol_for->sol, tab);
4026 int isl_basic_map_foreach_lexopt(__isl_keep isl_basic_map *bmap, int max,
4027 int (*fn)(__isl_take isl_basic_set *dom, __isl_take isl_mat *map,
4028 void *user),
4029 void *user)
4031 struct isl_sol_for *sol_for = NULL;
4033 bmap = isl_basic_map_copy(bmap);
4034 if (!bmap)
4035 return -1;
4037 bmap = isl_basic_map_detect_equalities(bmap);
4038 sol_for = sol_for_init(bmap, max, fn, user);
4040 if (isl_basic_map_fast_is_empty(bmap))
4041 /* nothing */;
4042 else {
4043 struct isl_tab *tab;
4044 struct isl_context *context = sol_for->sol.context;
4045 tab = tab_for_lexmin(bmap,
4046 context->op->peek_basic_set(context), 1, max);
4047 tab = context->op->detect_nonnegative_parameters(context, tab);
4048 sol_for_find_solutions(sol_for, tab);
4049 if (sol_for->sol.error)
4050 goto error;
4053 sol_free(&sol_for->sol);
4054 isl_basic_map_free(bmap);
4055 return 0;
4056 error:
4057 sol_free(&sol_for->sol);
4058 isl_basic_map_free(bmap);
4059 return -1;
4062 int isl_basic_map_foreach_lexmin(__isl_keep isl_basic_map *bmap,
4063 int (*fn)(__isl_take isl_basic_set *dom, __isl_take isl_mat *map,
4064 void *user),
4065 void *user)
4067 return isl_basic_map_foreach_lexopt(bmap, 0, fn, user);
4070 int isl_basic_map_foreach_lexmax(__isl_keep isl_basic_map *bmap,
4071 int (*fn)(__isl_take isl_basic_set *dom, __isl_take isl_mat *map,
4072 void *user),
4073 void *user)
4075 return isl_basic_map_foreach_lexopt(bmap, 1, fn, user);