isl_scheduler.c: compute_weights: use isl_basic_map_n_equality
[isl.git] / isl_tab_pip.c
blobf342a168f1815f528b05b9c26139be4b575c0a03
1 /*
2 * Copyright 2008-2009 Katholieke Universiteit Leuven
3 * Copyright 2010 INRIA Saclay
4 * Copyright 2016 Sven Verdoolaege
6 * Use of this software is governed by the MIT license
8 * Written by Sven Verdoolaege, K.U.Leuven, Departement
9 * Computerwetenschappen, Celestijnenlaan 200A, B-3001 Leuven, Belgium
10 * and INRIA Saclay - Ile-de-France, Parc Club Orsay Universite,
11 * ZAC des vignes, 4 rue Jacques Monod, 91893 Orsay, France
14 #include <isl_ctx_private.h>
15 #include "isl_map_private.h"
16 #include <isl_seq.h>
17 #include "isl_tab.h"
18 #include "isl_sample.h"
19 #include <isl_mat_private.h>
20 #include <isl_vec_private.h>
21 #include <isl_aff_private.h>
22 #include <isl_constraint_private.h>
23 #include <isl_options_private.h>
24 #include <isl_config.h>
26 #include <bset_to_bmap.c>
29 * The implementation of parametric integer linear programming in this file
30 * was inspired by the paper "Parametric Integer Programming" and the
31 * report "Solving systems of affine (in)equalities" by Paul Feautrier
32 * (and others).
34 * The strategy used for obtaining a feasible solution is different
35 * from the one used in isl_tab.c. In particular, in isl_tab.c,
36 * upon finding a constraint that is not yet satisfied, we pivot
37 * in a row that increases the constant term of the row holding the
38 * constraint, making sure the sample solution remains feasible
39 * for all the constraints it already satisfied.
40 * Here, we always pivot in the row holding the constraint,
41 * choosing a column that induces the lexicographically smallest
42 * increment to the sample solution.
44 * By starting out from a sample value that is lexicographically
45 * smaller than any integer point in the problem space, the first
46 * feasible integer sample point we find will also be the lexicographically
47 * smallest. If all variables can be assumed to be non-negative,
48 * then the initial sample value may be chosen equal to zero.
49 * However, we will not make this assumption. Instead, we apply
50 * the "big parameter" trick. Any variable x is then not directly
51 * used in the tableau, but instead it is represented by another
52 * variable x' = M + x, where M is an arbitrarily large (positive)
53 * value. x' is therefore always non-negative, whatever the value of x.
54 * Taking as initial sample value x' = 0 corresponds to x = -M,
55 * which is always smaller than any possible value of x.
57 * The big parameter trick is used in the main tableau and
58 * also in the context tableau if isl_context_lex is used.
59 * In this case, each tableaus has its own big parameter.
60 * Before doing any real work, we check if all the parameters
61 * happen to be non-negative. If so, we drop the column corresponding
62 * to M from the initial context tableau.
63 * If isl_context_gbr is used, then the big parameter trick is only
64 * used in the main tableau.
67 struct isl_context;
68 struct isl_context_op {
69 /* detect nonnegative parameters in context and mark them in tab */
70 struct isl_tab *(*detect_nonnegative_parameters)(
71 struct isl_context *context, struct isl_tab *tab);
72 /* return temporary reference to basic set representation of context */
73 struct isl_basic_set *(*peek_basic_set)(struct isl_context *context);
74 /* return temporary reference to tableau representation of context */
75 struct isl_tab *(*peek_tab)(struct isl_context *context);
76 /* add equality; check is 1 if eq may not be valid;
77 * update is 1 if we may want to call ineq_sign on context later.
79 void (*add_eq)(struct isl_context *context, isl_int *eq,
80 int check, int update);
81 /* add inequality; check is 1 if ineq may not be valid;
82 * update is 1 if we may want to call ineq_sign on context later.
84 void (*add_ineq)(struct isl_context *context, isl_int *ineq,
85 int check, int update);
86 /* check sign of ineq based on previous information.
87 * strict is 1 if saturation should be treated as a positive sign.
89 enum isl_tab_row_sign (*ineq_sign)(struct isl_context *context,
90 isl_int *ineq, int strict);
91 /* check if inequality maintains feasibility */
92 int (*test_ineq)(struct isl_context *context, isl_int *ineq);
93 /* return index of a div that corresponds to "div" */
94 int (*get_div)(struct isl_context *context, struct isl_tab *tab,
95 struct isl_vec *div);
96 /* insert div "div" to context at "pos" and return non-negativity */
97 isl_bool (*insert_div)(struct isl_context *context, int pos,
98 __isl_keep isl_vec *div);
99 int (*detect_equalities)(struct isl_context *context,
100 struct isl_tab *tab);
101 /* return row index of "best" split */
102 int (*best_split)(struct isl_context *context, struct isl_tab *tab);
103 /* check if context has already been determined to be empty */
104 int (*is_empty)(struct isl_context *context);
105 /* check if context is still usable */
106 int (*is_ok)(struct isl_context *context);
107 /* save a copy/snapshot of context */
108 void *(*save)(struct isl_context *context);
109 /* restore saved context */
110 void (*restore)(struct isl_context *context, void *);
111 /* discard saved context */
112 void (*discard)(void *);
113 /* invalidate context */
114 void (*invalidate)(struct isl_context *context);
115 /* free context */
116 __isl_null struct isl_context *(*free)(struct isl_context *context);
119 /* Shared parts of context representation.
121 * "n_unknown" is the number of final unknown integer divisions
122 * in the input domain.
124 struct isl_context {
125 struct isl_context_op *op;
126 int n_unknown;
129 struct isl_context_lex {
130 struct isl_context context;
131 struct isl_tab *tab;
134 /* A stack (linked list) of solutions of subtrees of the search space.
136 * "ma" describes the solution as a function of "dom".
137 * In particular, the domain space of "ma" is equal to the space of "dom".
139 * If "ma" is NULL, then there is no solution on "dom".
141 struct isl_partial_sol {
142 int level;
143 struct isl_basic_set *dom;
144 isl_multi_aff *ma;
146 struct isl_partial_sol *next;
149 struct isl_sol;
150 struct isl_sol_callback {
151 struct isl_tab_callback callback;
152 struct isl_sol *sol;
155 /* isl_sol is an interface for constructing a solution to
156 * a parametric integer linear programming problem.
157 * Every time the algorithm reaches a state where a solution
158 * can be read off from the tableau, the function "add" is called
159 * on the isl_sol passed to find_solutions_main. In a state where
160 * the tableau is empty, "add_empty" is called instead.
161 * "free" is called to free the implementation specific fields, if any.
163 * "error" is set if some error has occurred. This flag invalidates
164 * the remainder of the data structure.
165 * If "rational" is set, then a rational optimization is being performed.
166 * "level" is the current level in the tree with nodes for each
167 * split in the context.
168 * If "max" is set, then a maximization problem is being solved, rather than
169 * a minimization problem, which means that the variables in the
170 * tableau have value "M - x" rather than "M + x".
171 * "n_out" is the number of output dimensions in the input.
172 * "space" is the space in which the solution (and also the input) lives.
174 * The context tableau is owned by isl_sol and is updated incrementally.
176 * There are currently three implementations of this interface,
177 * isl_sol_map, which simply collects the solutions in an isl_map
178 * and (optionally) the parts of the context where there is no solution
179 * in an isl_set,
180 * isl_sol_pma, which collects an isl_pw_multi_aff instead, and
181 * isl_sol_for, which calls a user-defined function for each part of
182 * the solution.
184 struct isl_sol {
185 int error;
186 int rational;
187 int level;
188 int max;
189 int n_out;
190 isl_space *space;
191 struct isl_context *context;
192 struct isl_partial_sol *partial;
193 void (*add)(struct isl_sol *sol,
194 __isl_take isl_basic_set *dom, __isl_take isl_multi_aff *ma);
195 void (*add_empty)(struct isl_sol *sol, struct isl_basic_set *bset);
196 void (*free)(struct isl_sol *sol);
197 struct isl_sol_callback dec_level;
200 static void sol_free(struct isl_sol *sol)
202 struct isl_partial_sol *partial, *next;
203 if (!sol)
204 return;
205 for (partial = sol->partial; partial; partial = next) {
206 next = partial->next;
207 isl_basic_set_free(partial->dom);
208 isl_multi_aff_free(partial->ma);
209 free(partial);
211 isl_space_free(sol->space);
212 if (sol->context)
213 sol->context->op->free(sol->context);
214 sol->free(sol);
215 free(sol);
218 /* Push a partial solution represented by a domain and function "ma"
219 * onto the stack of partial solutions.
220 * If "ma" is NULL, then "dom" represents a part of the domain
221 * with no solution.
223 static void sol_push_sol(struct isl_sol *sol,
224 __isl_take isl_basic_set *dom, __isl_take isl_multi_aff *ma)
226 struct isl_partial_sol *partial;
228 if (sol->error || !dom)
229 goto error;
231 partial = isl_alloc_type(dom->ctx, struct isl_partial_sol);
232 if (!partial)
233 goto error;
235 partial->level = sol->level;
236 partial->dom = dom;
237 partial->ma = ma;
238 partial->next = sol->partial;
240 sol->partial = partial;
242 return;
243 error:
244 isl_basic_set_free(dom);
245 isl_multi_aff_free(ma);
246 sol->error = 1;
249 /* Check that the final columns of "M", starting at "first", are zero.
251 static isl_stat check_final_columns_are_zero(__isl_keep isl_mat *M,
252 unsigned first)
254 int i;
255 unsigned rows, cols, n;
257 if (!M)
258 return isl_stat_error;
259 rows = isl_mat_rows(M);
260 cols = isl_mat_cols(M);
261 n = cols - first;
262 for (i = 0; i < rows; ++i)
263 if (isl_seq_first_non_zero(M->row[i] + first, n) != -1)
264 isl_die(isl_mat_get_ctx(M), isl_error_internal,
265 "final columns should be zero",
266 return isl_stat_error);
267 return isl_stat_ok;
270 /* Set the affine expressions in "ma" according to the rows in "M", which
271 * are defined over the local space "ls".
272 * The matrix "M" may have extra (zero) columns beyond the number
273 * of variables in "ls".
275 static __isl_give isl_multi_aff *set_from_affine_matrix(
276 __isl_take isl_multi_aff *ma, __isl_take isl_local_space *ls,
277 __isl_take isl_mat *M)
279 int i, dim;
280 isl_aff *aff;
282 if (!ma || !ls || !M)
283 goto error;
285 dim = isl_local_space_dim(ls, isl_dim_all);
286 if (check_final_columns_are_zero(M, 1 + dim) < 0)
287 goto error;
288 for (i = 1; i < M->n_row; ++i) {
289 aff = isl_aff_alloc(isl_local_space_copy(ls));
290 if (aff) {
291 isl_int_set(aff->v->el[0], M->row[0][0]);
292 isl_seq_cpy(aff->v->el + 1, M->row[i], 1 + dim);
294 aff = isl_aff_normalize(aff);
295 ma = isl_multi_aff_set_aff(ma, i - 1, aff);
297 isl_local_space_free(ls);
298 isl_mat_free(M);
300 return ma;
301 error:
302 isl_local_space_free(ls);
303 isl_mat_free(M);
304 isl_multi_aff_free(ma);
305 return NULL;
308 /* Push a partial solution represented by a domain and mapping M
309 * onto the stack of partial solutions.
311 * The affine matrix "M" maps the dimensions of the context
312 * to the output variables. Convert it into an isl_multi_aff and
313 * then call sol_push_sol.
315 * Note that the description of the initial context may have involved
316 * existentially quantified variables, in which case they also appear
317 * in "dom". These need to be removed before creating the affine
318 * expression because an affine expression cannot be defined in terms
319 * of existentially quantified variables without a known representation.
320 * Since newly added integer divisions are inserted before these
321 * existentially quantified variables, they are still in the final
322 * positions and the corresponding final columns of "M" are zero
323 * because align_context_divs adds the existentially quantified
324 * variables of the context to the main tableau without any constraints and
325 * any equality constraints that are added later on can only serve
326 * to eliminate these existentially quantified variables.
328 static void sol_push_sol_mat(struct isl_sol *sol,
329 __isl_take isl_basic_set *dom, __isl_take isl_mat *M)
331 isl_local_space *ls;
332 isl_multi_aff *ma;
333 int n_div, n_known;
335 n_div = isl_basic_set_dim(dom, isl_dim_div);
336 n_known = n_div - sol->context->n_unknown;
338 ma = isl_multi_aff_alloc(isl_space_copy(sol->space));
339 ls = isl_basic_set_get_local_space(dom);
340 ls = isl_local_space_drop_dims(ls, isl_dim_div,
341 n_known, n_div - n_known);
342 ma = set_from_affine_matrix(ma, ls, M);
344 if (!ma)
345 dom = isl_basic_set_free(dom);
346 sol_push_sol(sol, dom, ma);
349 /* Pop one partial solution from the partial solution stack and
350 * pass it on to sol->add or sol->add_empty.
352 static void sol_pop_one(struct isl_sol *sol)
354 struct isl_partial_sol *partial;
356 partial = sol->partial;
357 sol->partial = partial->next;
359 if (partial->ma)
360 sol->add(sol, partial->dom, partial->ma);
361 else
362 sol->add_empty(sol, partial->dom);
363 free(partial);
366 /* Return a fresh copy of the domain represented by the context tableau.
368 static struct isl_basic_set *sol_domain(struct isl_sol *sol)
370 struct isl_basic_set *bset;
372 if (sol->error)
373 return NULL;
375 bset = isl_basic_set_dup(sol->context->op->peek_basic_set(sol->context));
376 bset = isl_basic_set_update_from_tab(bset,
377 sol->context->op->peek_tab(sol->context));
379 return bset;
382 /* Check whether two partial solutions have the same affine expressions.
384 static isl_bool same_solution(struct isl_partial_sol *s1,
385 struct isl_partial_sol *s2)
387 if (!s1->ma != !s2->ma)
388 return isl_bool_false;
389 if (!s1->ma)
390 return isl_bool_true;
392 return isl_multi_aff_plain_is_equal(s1->ma, s2->ma);
395 /* Swap the initial two partial solutions in "sol".
397 * That is, go from
399 * sol->partial = p1; p1->next = p2; p2->next = p3
401 * to
403 * sol->partial = p2; p2->next = p1; p1->next = p3
405 static void swap_initial(struct isl_sol *sol)
407 struct isl_partial_sol *partial;
409 partial = sol->partial;
410 sol->partial = partial->next;
411 partial->next = partial->next->next;
412 sol->partial->next = partial;
415 /* Combine the initial two partial solution of "sol" into
416 * a partial solution with the current context domain of "sol" and
417 * the function description of the second partial solution in the list.
418 * The level of the new partial solution is set to the current level.
420 * That is, the first two partial solutions (D1,M1) and (D2,M2) are
421 * replaced by (D,M2), where D is the domain of "sol", which is assumed
422 * to be the union of D1 and D2, while M1 is assumed to be equal to M2
423 * (at least on D1).
425 static isl_stat combine_initial_into_second(struct isl_sol *sol)
427 struct isl_partial_sol *partial;
428 isl_basic_set *bset;
430 partial = sol->partial;
432 bset = sol_domain(sol);
433 isl_basic_set_free(partial->next->dom);
434 partial->next->dom = bset;
435 partial->next->level = sol->level;
437 if (!bset)
438 return isl_stat_error;
440 sol->partial = partial->next;
441 isl_basic_set_free(partial->dom);
442 isl_multi_aff_free(partial->ma);
443 free(partial);
445 return isl_stat_ok;
448 /* Are "ma1" and "ma2" equal to each other on "dom"?
450 * Combine "ma1" and "ma2" with "dom" and check if the results are the same.
451 * "dom" may have existentially quantified variables. Eliminate them first
452 * as otherwise they would have to be eliminated twice, in a more complicated
453 * context.
455 static isl_bool equal_on_domain(__isl_keep isl_multi_aff *ma1,
456 __isl_keep isl_multi_aff *ma2, __isl_keep isl_basic_set *dom)
458 isl_set *set;
459 isl_pw_multi_aff *pma1, *pma2;
460 isl_bool equal;
462 set = isl_basic_set_compute_divs(isl_basic_set_copy(dom));
463 pma1 = isl_pw_multi_aff_alloc(isl_set_copy(set),
464 isl_multi_aff_copy(ma1));
465 pma2 = isl_pw_multi_aff_alloc(set, isl_multi_aff_copy(ma2));
466 equal = isl_pw_multi_aff_is_equal(pma1, pma2);
467 isl_pw_multi_aff_free(pma1);
468 isl_pw_multi_aff_free(pma2);
470 return equal;
473 /* The initial two partial solutions of "sol" are known to be at
474 * the same level.
475 * If they represent the same solution (on different parts of the domain),
476 * then combine them into a single solution at the current level.
477 * Otherwise, pop them both.
479 * Even if the two partial solution are not obviously the same,
480 * one may still be a simplification of the other over its own domain.
481 * Also check if the two sets of affine functions are equal when
482 * restricted to one of the domains. If so, combine the two
483 * using the set of affine functions on the other domain.
484 * That is, for two partial solutions (D1,M1) and (D2,M2),
485 * if M1 = M2 on D1, then the pair of partial solutions can
486 * be replaced by (D1+D2,M2) and similarly when M1 = M2 on D2.
488 static isl_stat combine_initial_if_equal(struct isl_sol *sol)
490 struct isl_partial_sol *partial;
491 isl_bool same;
493 partial = sol->partial;
495 same = same_solution(partial, partial->next);
496 if (same < 0)
497 return isl_stat_error;
498 if (same)
499 return combine_initial_into_second(sol);
500 if (partial->ma && partial->next->ma) {
501 same = equal_on_domain(partial->ma, partial->next->ma,
502 partial->dom);
503 if (same < 0)
504 return isl_stat_error;
505 if (same)
506 return combine_initial_into_second(sol);
507 same = equal_on_domain(partial->ma, partial->next->ma,
508 partial->next->dom);
509 if (same) {
510 swap_initial(sol);
511 return combine_initial_into_second(sol);
515 sol_pop_one(sol);
516 sol_pop_one(sol);
518 return isl_stat_ok;
521 /* Pop all solutions from the partial solution stack that were pushed onto
522 * the stack at levels that are deeper than the current level.
523 * If the two topmost elements on the stack have the same level
524 * and represent the same solution, then their domains are combined.
525 * This combined domain is the same as the current context domain
526 * as sol_pop is called each time we move back to a higher level.
527 * If the outer level (0) has been reached, then all partial solutions
528 * at the current level are also popped off.
530 static void sol_pop(struct isl_sol *sol)
532 struct isl_partial_sol *partial;
534 if (sol->error)
535 return;
537 partial = sol->partial;
538 if (!partial)
539 return;
541 if (partial->level == 0 && sol->level == 0) {
542 for (partial = sol->partial; partial; partial = sol->partial)
543 sol_pop_one(sol);
544 return;
547 if (partial->level <= sol->level)
548 return;
550 if (partial->next && partial->next->level == partial->level) {
551 if (combine_initial_if_equal(sol) < 0)
552 goto error;
553 } else
554 sol_pop_one(sol);
556 if (sol->level == 0) {
557 for (partial = sol->partial; partial; partial = sol->partial)
558 sol_pop_one(sol);
559 return;
562 if (0)
563 error: sol->error = 1;
566 static void sol_dec_level(struct isl_sol *sol)
568 if (sol->error)
569 return;
571 sol->level--;
573 sol_pop(sol);
576 static isl_stat sol_dec_level_wrap(struct isl_tab_callback *cb)
578 struct isl_sol_callback *callback = (struct isl_sol_callback *)cb;
580 sol_dec_level(callback->sol);
582 return callback->sol->error ? isl_stat_error : isl_stat_ok;
585 /* Move down to next level and push callback onto context tableau
586 * to decrease the level again when it gets rolled back across
587 * the current state. That is, dec_level will be called with
588 * the context tableau in the same state as it is when inc_level
589 * is called.
591 static void sol_inc_level(struct isl_sol *sol)
593 struct isl_tab *tab;
595 if (sol->error)
596 return;
598 sol->level++;
599 tab = sol->context->op->peek_tab(sol->context);
600 if (isl_tab_push_callback(tab, &sol->dec_level.callback) < 0)
601 sol->error = 1;
604 static void scale_rows(struct isl_mat *mat, isl_int m, int n_row)
606 int i;
608 if (isl_int_is_one(m))
609 return;
611 for (i = 0; i < n_row; ++i)
612 isl_seq_scale(mat->row[i], mat->row[i], m, mat->n_col);
615 /* Add the solution identified by the tableau and the context tableau.
617 * The layout of the variables is as follows.
618 * tab->n_var is equal to the total number of variables in the input
619 * map (including divs that were copied from the context)
620 * + the number of extra divs constructed
621 * Of these, the first tab->n_param and the last tab->n_div variables
622 * correspond to the variables in the context, i.e.,
623 * tab->n_param + tab->n_div = context_tab->n_var
624 * tab->n_param is equal to the number of parameters and input
625 * dimensions in the input map
626 * tab->n_div is equal to the number of divs in the context
628 * If there is no solution, then call add_empty with a basic set
629 * that corresponds to the context tableau. (If add_empty is NULL,
630 * then do nothing).
632 * If there is a solution, then first construct a matrix that maps
633 * all dimensions of the context to the output variables, i.e.,
634 * the output dimensions in the input map.
635 * The divs in the input map (if any) that do not correspond to any
636 * div in the context do not appear in the solution.
637 * The algorithm will make sure that they have an integer value,
638 * but these values themselves are of no interest.
639 * We have to be careful not to drop or rearrange any divs in the
640 * context because that would change the meaning of the matrix.
642 * To extract the value of the output variables, it should be noted
643 * that we always use a big parameter M in the main tableau and so
644 * the variable stored in this tableau is not an output variable x itself, but
645 * x' = M + x (in case of minimization)
646 * or
647 * x' = M - x (in case of maximization)
648 * If x' appears in a column, then its optimal value is zero,
649 * which means that the optimal value of x is an unbounded number
650 * (-M for minimization and M for maximization).
651 * We currently assume that the output dimensions in the original map
652 * are bounded, so this cannot occur.
653 * Similarly, when x' appears in a row, then the coefficient of M in that
654 * row is necessarily 1.
655 * If the row in the tableau represents
656 * d x' = c + d M + e(y)
657 * then, in case of minimization, the corresponding row in the matrix
658 * will be
659 * a c + a e(y)
660 * with a d = m, the (updated) common denominator of the matrix.
661 * In case of maximization, the row will be
662 * -a c - a e(y)
664 static void sol_add(struct isl_sol *sol, struct isl_tab *tab)
666 struct isl_basic_set *bset = NULL;
667 struct isl_mat *mat = NULL;
668 unsigned off;
669 int row;
670 isl_int m;
672 if (sol->error || !tab)
673 goto error;
675 if (tab->empty && !sol->add_empty)
676 return;
677 if (sol->context->op->is_empty(sol->context))
678 return;
680 bset = sol_domain(sol);
682 if (tab->empty) {
683 sol_push_sol(sol, bset, NULL);
684 return;
687 off = 2 + tab->M;
689 mat = isl_mat_alloc(tab->mat->ctx, 1 + sol->n_out,
690 1 + tab->n_param + tab->n_div);
691 if (!mat)
692 goto error;
694 isl_int_init(m);
696 isl_seq_clr(mat->row[0] + 1, mat->n_col - 1);
697 isl_int_set_si(mat->row[0][0], 1);
698 for (row = 0; row < sol->n_out; ++row) {
699 int i = tab->n_param + row;
700 int r, j;
702 isl_seq_clr(mat->row[1 + row], mat->n_col);
703 if (!tab->var[i].is_row) {
704 if (tab->M)
705 isl_die(mat->ctx, isl_error_invalid,
706 "unbounded optimum", goto error2);
707 continue;
710 r = tab->var[i].index;
711 if (tab->M &&
712 isl_int_ne(tab->mat->row[r][2], tab->mat->row[r][0]))
713 isl_die(mat->ctx, isl_error_invalid,
714 "unbounded optimum", goto error2);
715 isl_int_gcd(m, mat->row[0][0], tab->mat->row[r][0]);
716 isl_int_divexact(m, tab->mat->row[r][0], m);
717 scale_rows(mat, m, 1 + row);
718 isl_int_divexact(m, mat->row[0][0], tab->mat->row[r][0]);
719 isl_int_mul(mat->row[1 + row][0], m, tab->mat->row[r][1]);
720 for (j = 0; j < tab->n_param; ++j) {
721 int col;
722 if (tab->var[j].is_row)
723 continue;
724 col = tab->var[j].index;
725 isl_int_mul(mat->row[1 + row][1 + j], m,
726 tab->mat->row[r][off + col]);
728 for (j = 0; j < tab->n_div; ++j) {
729 int col;
730 if (tab->var[tab->n_var - tab->n_div+j].is_row)
731 continue;
732 col = tab->var[tab->n_var - tab->n_div+j].index;
733 isl_int_mul(mat->row[1 + row][1 + tab->n_param + j], m,
734 tab->mat->row[r][off + col]);
736 if (sol->max)
737 isl_seq_neg(mat->row[1 + row], mat->row[1 + row],
738 mat->n_col);
741 isl_int_clear(m);
743 sol_push_sol_mat(sol, bset, mat);
744 return;
745 error2:
746 isl_int_clear(m);
747 error:
748 isl_basic_set_free(bset);
749 isl_mat_free(mat);
750 sol->error = 1;
753 struct isl_sol_map {
754 struct isl_sol sol;
755 struct isl_map *map;
756 struct isl_set *empty;
759 static void sol_map_free(struct isl_sol *sol)
761 struct isl_sol_map *sol_map = (struct isl_sol_map *) sol;
762 isl_map_free(sol_map->map);
763 isl_set_free(sol_map->empty);
766 /* This function is called for parts of the context where there is
767 * no solution, with "bset" corresponding to the context tableau.
768 * Simply add the basic set to the set "empty".
770 static void sol_map_add_empty(struct isl_sol_map *sol,
771 struct isl_basic_set *bset)
773 if (!bset || !sol->empty)
774 goto error;
776 sol->empty = isl_set_grow(sol->empty, 1);
777 bset = isl_basic_set_simplify(bset);
778 bset = isl_basic_set_finalize(bset);
779 sol->empty = isl_set_add_basic_set(sol->empty, isl_basic_set_copy(bset));
780 if (!sol->empty)
781 goto error;
782 isl_basic_set_free(bset);
783 return;
784 error:
785 isl_basic_set_free(bset);
786 sol->sol.error = 1;
789 static void sol_map_add_empty_wrap(struct isl_sol *sol,
790 struct isl_basic_set *bset)
792 sol_map_add_empty((struct isl_sol_map *)sol, bset);
795 /* Given a basic set "dom" that represents the context and a tuple of
796 * affine expressions "ma" defined over this domain, construct a basic map
797 * that expresses this function on the domain.
799 static void sol_map_add(struct isl_sol_map *sol,
800 __isl_take isl_basic_set *dom, __isl_take isl_multi_aff *ma)
802 isl_basic_map *bmap;
804 if (sol->sol.error || !dom || !ma)
805 goto error;
807 bmap = isl_basic_map_from_multi_aff2(ma, sol->sol.rational);
808 bmap = isl_basic_map_intersect_domain(bmap, dom);
809 sol->map = isl_map_grow(sol->map, 1);
810 sol->map = isl_map_add_basic_map(sol->map, bmap);
811 if (!sol->map)
812 sol->sol.error = 1;
813 return;
814 error:
815 isl_basic_set_free(dom);
816 isl_multi_aff_free(ma);
817 sol->sol.error = 1;
820 static void sol_map_add_wrap(struct isl_sol *sol,
821 __isl_take isl_basic_set *dom, __isl_take isl_multi_aff *ma)
823 sol_map_add((struct isl_sol_map *)sol, dom, ma);
827 /* Store the "parametric constant" of row "row" of tableau "tab" in "line",
828 * i.e., the constant term and the coefficients of all variables that
829 * appear in the context tableau.
830 * Note that the coefficient of the big parameter M is NOT copied.
831 * The context tableau may not have a big parameter and even when it
832 * does, it is a different big parameter.
834 static void get_row_parameter_line(struct isl_tab *tab, int row, isl_int *line)
836 int i;
837 unsigned off = 2 + tab->M;
839 isl_int_set(line[0], tab->mat->row[row][1]);
840 for (i = 0; i < tab->n_param; ++i) {
841 if (tab->var[i].is_row)
842 isl_int_set_si(line[1 + i], 0);
843 else {
844 int col = tab->var[i].index;
845 isl_int_set(line[1 + i], tab->mat->row[row][off + col]);
848 for (i = 0; i < tab->n_div; ++i) {
849 if (tab->var[tab->n_var - tab->n_div + i].is_row)
850 isl_int_set_si(line[1 + tab->n_param + i], 0);
851 else {
852 int col = tab->var[tab->n_var - tab->n_div + i].index;
853 isl_int_set(line[1 + tab->n_param + i],
854 tab->mat->row[row][off + col]);
859 /* Check if rows "row1" and "row2" have identical "parametric constants",
860 * as explained above.
861 * In this case, we also insist that the coefficients of the big parameter
862 * be the same as the values of the constants will only be the same
863 * if these coefficients are also the same.
865 static int identical_parameter_line(struct isl_tab *tab, int row1, int row2)
867 int i;
868 unsigned off = 2 + tab->M;
870 if (isl_int_ne(tab->mat->row[row1][1], tab->mat->row[row2][1]))
871 return 0;
873 if (tab->M && isl_int_ne(tab->mat->row[row1][2],
874 tab->mat->row[row2][2]))
875 return 0;
877 for (i = 0; i < tab->n_param + tab->n_div; ++i) {
878 int pos = i < tab->n_param ? i :
879 tab->n_var - tab->n_div + i - tab->n_param;
880 int col;
882 if (tab->var[pos].is_row)
883 continue;
884 col = tab->var[pos].index;
885 if (isl_int_ne(tab->mat->row[row1][off + col],
886 tab->mat->row[row2][off + col]))
887 return 0;
889 return 1;
892 /* Return an inequality that expresses that the "parametric constant"
893 * should be non-negative.
894 * This function is only called when the coefficient of the big parameter
895 * is equal to zero.
897 static struct isl_vec *get_row_parameter_ineq(struct isl_tab *tab, int row)
899 struct isl_vec *ineq;
901 ineq = isl_vec_alloc(tab->mat->ctx, 1 + tab->n_param + tab->n_div);
902 if (!ineq)
903 return NULL;
905 get_row_parameter_line(tab, row, ineq->el);
906 if (ineq)
907 ineq = isl_vec_normalize(ineq);
909 return ineq;
912 /* Normalize a div expression of the form
914 * [(g*f(x) + c)/(g * m)]
916 * with c the constant term and f(x) the remaining coefficients, to
918 * [(f(x) + [c/g])/m]
920 static void normalize_div(__isl_keep isl_vec *div)
922 isl_ctx *ctx = isl_vec_get_ctx(div);
923 int len = div->size - 2;
925 isl_seq_gcd(div->el + 2, len, &ctx->normalize_gcd);
926 isl_int_gcd(ctx->normalize_gcd, ctx->normalize_gcd, div->el[0]);
928 if (isl_int_is_one(ctx->normalize_gcd))
929 return;
931 isl_int_divexact(div->el[0], div->el[0], ctx->normalize_gcd);
932 isl_int_fdiv_q(div->el[1], div->el[1], ctx->normalize_gcd);
933 isl_seq_scale_down(div->el + 2, div->el + 2, ctx->normalize_gcd, len);
936 /* Return an integer division for use in a parametric cut based
937 * on the given row.
938 * In particular, let the parametric constant of the row be
940 * \sum_i a_i y_i
942 * where y_0 = 1, but none of the y_i corresponds to the big parameter M.
943 * The div returned is equal to
945 * floor(\sum_i {-a_i} y_i) = floor((\sum_i (-a_i mod d) y_i)/d)
947 static struct isl_vec *get_row_parameter_div(struct isl_tab *tab, int row)
949 struct isl_vec *div;
951 div = isl_vec_alloc(tab->mat->ctx, 1 + 1 + tab->n_param + tab->n_div);
952 if (!div)
953 return NULL;
955 isl_int_set(div->el[0], tab->mat->row[row][0]);
956 get_row_parameter_line(tab, row, div->el + 1);
957 isl_seq_neg(div->el + 1, div->el + 1, div->size - 1);
958 normalize_div(div);
959 isl_seq_fdiv_r(div->el + 1, div->el + 1, div->el[0], div->size - 1);
961 return div;
964 /* Return an integer division for use in transferring an integrality constraint
965 * to the context.
966 * In particular, let the parametric constant of the row be
968 * \sum_i a_i y_i
970 * where y_0 = 1, but none of the y_i corresponds to the big parameter M.
971 * The the returned div is equal to
973 * floor(\sum_i {a_i} y_i) = floor((\sum_i (a_i mod d) y_i)/d)
975 static struct isl_vec *get_row_split_div(struct isl_tab *tab, int row)
977 struct isl_vec *div;
979 div = isl_vec_alloc(tab->mat->ctx, 1 + 1 + tab->n_param + tab->n_div);
980 if (!div)
981 return NULL;
983 isl_int_set(div->el[0], tab->mat->row[row][0]);
984 get_row_parameter_line(tab, row, div->el + 1);
985 normalize_div(div);
986 isl_seq_fdiv_r(div->el + 1, div->el + 1, div->el[0], div->size - 1);
988 return div;
991 /* Construct and return an inequality that expresses an upper bound
992 * on the given div.
993 * In particular, if the div is given by
995 * d = floor(e/m)
997 * then the inequality expresses
999 * m d <= e
1001 static __isl_give isl_vec *ineq_for_div(__isl_keep isl_basic_set *bset,
1002 unsigned div)
1004 unsigned total;
1005 unsigned div_pos;
1006 struct isl_vec *ineq;
1008 if (!bset)
1009 return NULL;
1011 total = isl_basic_set_total_dim(bset);
1012 div_pos = 1 + total - bset->n_div + div;
1014 ineq = isl_vec_alloc(bset->ctx, 1 + total);
1015 if (!ineq)
1016 return NULL;
1018 isl_seq_cpy(ineq->el, bset->div[div] + 1, 1 + total);
1019 isl_int_neg(ineq->el[div_pos], bset->div[div][0]);
1020 return ineq;
1023 /* Given a row in the tableau and a div that was created
1024 * using get_row_split_div and that has been constrained to equality, i.e.,
1026 * d = floor(\sum_i {a_i} y_i) = \sum_i {a_i} y_i
1028 * replace the expression "\sum_i {a_i} y_i" in the row by d,
1029 * i.e., we subtract "\sum_i {a_i} y_i" and add 1 d.
1030 * The coefficients of the non-parameters in the tableau have been
1031 * verified to be integral. We can therefore simply replace coefficient b
1032 * by floor(b). For the coefficients of the parameters we have
1033 * floor(a_i) = a_i - {a_i}, while for the other coefficients, we have
1034 * floor(b) = b.
1036 static struct isl_tab *set_row_cst_to_div(struct isl_tab *tab, int row, int div)
1038 isl_seq_fdiv_q(tab->mat->row[row] + 1, tab->mat->row[row] + 1,
1039 tab->mat->row[row][0], 1 + tab->M + tab->n_col);
1041 isl_int_set_si(tab->mat->row[row][0], 1);
1043 if (tab->var[tab->n_var - tab->n_div + div].is_row) {
1044 int drow = tab->var[tab->n_var - tab->n_div + div].index;
1046 isl_assert(tab->mat->ctx,
1047 isl_int_is_one(tab->mat->row[drow][0]), goto error);
1048 isl_seq_combine(tab->mat->row[row] + 1,
1049 tab->mat->ctx->one, tab->mat->row[row] + 1,
1050 tab->mat->ctx->one, tab->mat->row[drow] + 1,
1051 1 + tab->M + tab->n_col);
1052 } else {
1053 int dcol = tab->var[tab->n_var - tab->n_div + div].index;
1055 isl_int_add_ui(tab->mat->row[row][2 + tab->M + dcol],
1056 tab->mat->row[row][2 + tab->M + dcol], 1);
1059 return tab;
1060 error:
1061 isl_tab_free(tab);
1062 return NULL;
1065 /* Check if the (parametric) constant of the given row is obviously
1066 * negative, meaning that we don't need to consult the context tableau.
1067 * If there is a big parameter and its coefficient is non-zero,
1068 * then this coefficient determines the outcome.
1069 * Otherwise, we check whether the constant is negative and
1070 * all non-zero coefficients of parameters are negative and
1071 * belong to non-negative parameters.
1073 static int is_obviously_neg(struct isl_tab *tab, int row)
1075 int i;
1076 int col;
1077 unsigned off = 2 + tab->M;
1079 if (tab->M) {
1080 if (isl_int_is_pos(tab->mat->row[row][2]))
1081 return 0;
1082 if (isl_int_is_neg(tab->mat->row[row][2]))
1083 return 1;
1086 if (isl_int_is_nonneg(tab->mat->row[row][1]))
1087 return 0;
1088 for (i = 0; i < tab->n_param; ++i) {
1089 /* Eliminated parameter */
1090 if (tab->var[i].is_row)
1091 continue;
1092 col = tab->var[i].index;
1093 if (isl_int_is_zero(tab->mat->row[row][off + col]))
1094 continue;
1095 if (!tab->var[i].is_nonneg)
1096 return 0;
1097 if (isl_int_is_pos(tab->mat->row[row][off + col]))
1098 return 0;
1100 for (i = 0; i < tab->n_div; ++i) {
1101 if (tab->var[tab->n_var - tab->n_div + i].is_row)
1102 continue;
1103 col = tab->var[tab->n_var - tab->n_div + i].index;
1104 if (isl_int_is_zero(tab->mat->row[row][off + col]))
1105 continue;
1106 if (!tab->var[tab->n_var - tab->n_div + i].is_nonneg)
1107 return 0;
1108 if (isl_int_is_pos(tab->mat->row[row][off + col]))
1109 return 0;
1111 return 1;
1114 /* Check if the (parametric) constant of the given row is obviously
1115 * non-negative, meaning that we don't need to consult the context tableau.
1116 * If there is a big parameter and its coefficient is non-zero,
1117 * then this coefficient determines the outcome.
1118 * Otherwise, we check whether the constant is non-negative and
1119 * all non-zero coefficients of parameters are positive and
1120 * belong to non-negative parameters.
1122 static int is_obviously_nonneg(struct isl_tab *tab, int row)
1124 int i;
1125 int col;
1126 unsigned off = 2 + tab->M;
1128 if (tab->M) {
1129 if (isl_int_is_pos(tab->mat->row[row][2]))
1130 return 1;
1131 if (isl_int_is_neg(tab->mat->row[row][2]))
1132 return 0;
1135 if (isl_int_is_neg(tab->mat->row[row][1]))
1136 return 0;
1137 for (i = 0; i < tab->n_param; ++i) {
1138 /* Eliminated parameter */
1139 if (tab->var[i].is_row)
1140 continue;
1141 col = tab->var[i].index;
1142 if (isl_int_is_zero(tab->mat->row[row][off + col]))
1143 continue;
1144 if (!tab->var[i].is_nonneg)
1145 return 0;
1146 if (isl_int_is_neg(tab->mat->row[row][off + col]))
1147 return 0;
1149 for (i = 0; i < tab->n_div; ++i) {
1150 if (tab->var[tab->n_var - tab->n_div + i].is_row)
1151 continue;
1152 col = tab->var[tab->n_var - tab->n_div + i].index;
1153 if (isl_int_is_zero(tab->mat->row[row][off + col]))
1154 continue;
1155 if (!tab->var[tab->n_var - tab->n_div + i].is_nonneg)
1156 return 0;
1157 if (isl_int_is_neg(tab->mat->row[row][off + col]))
1158 return 0;
1160 return 1;
1163 /* Given a row r and two columns, return the column that would
1164 * lead to the lexicographically smallest increment in the sample
1165 * solution when leaving the basis in favor of the row.
1166 * Pivoting with column c will increment the sample value by a non-negative
1167 * constant times a_{V,c}/a_{r,c}, with a_{V,c} the elements of column c
1168 * corresponding to the non-parametric variables.
1169 * If variable v appears in a column c_v, the a_{v,c} = 1 iff c = c_v,
1170 * with all other entries in this virtual row equal to zero.
1171 * If variable v appears in a row, then a_{v,c} is the element in column c
1172 * of that row.
1174 * Let v be the first variable with a_{v,c1}/a_{r,c1} != a_{v,c2}/a_{r,c2}.
1175 * Then if a_{v,c1}/a_{r,c1} < a_{v,c2}/a_{r,c2}, i.e.,
1176 * a_{v,c2} a_{r,c1} - a_{v,c1} a_{r,c2} > 0, c1 results in the minimal
1177 * increment. Otherwise, it's c2.
1179 static int lexmin_col_pair(struct isl_tab *tab,
1180 int row, int col1, int col2, isl_int tmp)
1182 int i;
1183 isl_int *tr;
1185 tr = tab->mat->row[row] + 2 + tab->M;
1187 for (i = tab->n_param; i < tab->n_var - tab->n_div; ++i) {
1188 int s1, s2;
1189 isl_int *r;
1191 if (!tab->var[i].is_row) {
1192 if (tab->var[i].index == col1)
1193 return col2;
1194 if (tab->var[i].index == col2)
1195 return col1;
1196 continue;
1199 if (tab->var[i].index == row)
1200 continue;
1202 r = tab->mat->row[tab->var[i].index] + 2 + tab->M;
1203 s1 = isl_int_sgn(r[col1]);
1204 s2 = isl_int_sgn(r[col2]);
1205 if (s1 == 0 && s2 == 0)
1206 continue;
1207 if (s1 < s2)
1208 return col1;
1209 if (s2 < s1)
1210 return col2;
1212 isl_int_mul(tmp, r[col2], tr[col1]);
1213 isl_int_submul(tmp, r[col1], tr[col2]);
1214 if (isl_int_is_pos(tmp))
1215 return col1;
1216 if (isl_int_is_neg(tmp))
1217 return col2;
1219 return -1;
1222 /* Given a row in the tableau, find and return the column that would
1223 * result in the lexicographically smallest, but positive, increment
1224 * in the sample point.
1225 * If there is no such column, then return tab->n_col.
1226 * If anything goes wrong, return -1.
1228 static int lexmin_pivot_col(struct isl_tab *tab, int row)
1230 int j;
1231 int col = tab->n_col;
1232 isl_int *tr;
1233 isl_int tmp;
1235 tr = tab->mat->row[row] + 2 + tab->M;
1237 isl_int_init(tmp);
1239 for (j = tab->n_dead; j < tab->n_col; ++j) {
1240 if (tab->col_var[j] >= 0 &&
1241 (tab->col_var[j] < tab->n_param ||
1242 tab->col_var[j] >= tab->n_var - tab->n_div))
1243 continue;
1245 if (!isl_int_is_pos(tr[j]))
1246 continue;
1248 if (col == tab->n_col)
1249 col = j;
1250 else
1251 col = lexmin_col_pair(tab, row, col, j, tmp);
1252 isl_assert(tab->mat->ctx, col >= 0, goto error);
1255 isl_int_clear(tmp);
1256 return col;
1257 error:
1258 isl_int_clear(tmp);
1259 return -1;
1262 /* Return the first known violated constraint, i.e., a non-negative
1263 * constraint that currently has an either obviously negative value
1264 * or a previously determined to be negative value.
1266 * If any constraint has a negative coefficient for the big parameter,
1267 * if any, then we return one of these first.
1269 static int first_neg(struct isl_tab *tab)
1271 int row;
1273 if (tab->M)
1274 for (row = tab->n_redundant; row < tab->n_row; ++row) {
1275 if (!isl_tab_var_from_row(tab, row)->is_nonneg)
1276 continue;
1277 if (!isl_int_is_neg(tab->mat->row[row][2]))
1278 continue;
1279 if (tab->row_sign)
1280 tab->row_sign[row] = isl_tab_row_neg;
1281 return row;
1283 for (row = tab->n_redundant; row < tab->n_row; ++row) {
1284 if (!isl_tab_var_from_row(tab, row)->is_nonneg)
1285 continue;
1286 if (tab->row_sign) {
1287 if (tab->row_sign[row] == 0 &&
1288 is_obviously_neg(tab, row))
1289 tab->row_sign[row] = isl_tab_row_neg;
1290 if (tab->row_sign[row] != isl_tab_row_neg)
1291 continue;
1292 } else if (!is_obviously_neg(tab, row))
1293 continue;
1294 return row;
1296 return -1;
1299 /* Check whether the invariant that all columns are lexico-positive
1300 * is satisfied. This function is not called from the current code
1301 * but is useful during debugging.
1303 static void check_lexpos(struct isl_tab *tab) __attribute__ ((unused));
1304 static void check_lexpos(struct isl_tab *tab)
1306 unsigned off = 2 + tab->M;
1307 int col;
1308 int var;
1309 int row;
1311 for (col = tab->n_dead; col < tab->n_col; ++col) {
1312 if (tab->col_var[col] >= 0 &&
1313 (tab->col_var[col] < tab->n_param ||
1314 tab->col_var[col] >= tab->n_var - tab->n_div))
1315 continue;
1316 for (var = tab->n_param; var < tab->n_var - tab->n_div; ++var) {
1317 if (!tab->var[var].is_row) {
1318 if (tab->var[var].index == col)
1319 break;
1320 else
1321 continue;
1323 row = tab->var[var].index;
1324 if (isl_int_is_zero(tab->mat->row[row][off + col]))
1325 continue;
1326 if (isl_int_is_pos(tab->mat->row[row][off + col]))
1327 break;
1328 fprintf(stderr, "lexneg column %d (row %d)\n",
1329 col, row);
1331 if (var >= tab->n_var - tab->n_div)
1332 fprintf(stderr, "zero column %d\n", col);
1336 /* Report to the caller that the given constraint is part of an encountered
1337 * conflict.
1339 static int report_conflicting_constraint(struct isl_tab *tab, int con)
1341 return tab->conflict(con, tab->conflict_user);
1344 /* Given a conflicting row in the tableau, report all constraints
1345 * involved in the row to the caller. That is, the row itself
1346 * (if it represents a constraint) and all constraint columns with
1347 * non-zero (and therefore negative) coefficients.
1349 static int report_conflict(struct isl_tab *tab, int row)
1351 int j;
1352 isl_int *tr;
1354 if (!tab->conflict)
1355 return 0;
1357 if (tab->row_var[row] < 0 &&
1358 report_conflicting_constraint(tab, ~tab->row_var[row]) < 0)
1359 return -1;
1361 tr = tab->mat->row[row] + 2 + tab->M;
1363 for (j = tab->n_dead; j < tab->n_col; ++j) {
1364 if (tab->col_var[j] >= 0 &&
1365 (tab->col_var[j] < tab->n_param ||
1366 tab->col_var[j] >= tab->n_var - tab->n_div))
1367 continue;
1369 if (!isl_int_is_neg(tr[j]))
1370 continue;
1372 if (tab->col_var[j] < 0 &&
1373 report_conflicting_constraint(tab, ~tab->col_var[j]) < 0)
1374 return -1;
1377 return 0;
1380 /* Resolve all known or obviously violated constraints through pivoting.
1381 * In particular, as long as we can find any violated constraint, we
1382 * look for a pivoting column that would result in the lexicographically
1383 * smallest increment in the sample point. If there is no such column
1384 * then the tableau is infeasible.
1386 static int restore_lexmin(struct isl_tab *tab) WARN_UNUSED;
1387 static int restore_lexmin(struct isl_tab *tab)
1389 int row, col;
1391 if (!tab)
1392 return -1;
1393 if (tab->empty)
1394 return 0;
1395 while ((row = first_neg(tab)) != -1) {
1396 col = lexmin_pivot_col(tab, row);
1397 if (col >= tab->n_col) {
1398 if (report_conflict(tab, row) < 0)
1399 return -1;
1400 if (isl_tab_mark_empty(tab) < 0)
1401 return -1;
1402 return 0;
1404 if (col < 0)
1405 return -1;
1406 if (isl_tab_pivot(tab, row, col) < 0)
1407 return -1;
1409 return 0;
1412 /* Given a row that represents an equality, look for an appropriate
1413 * pivoting column.
1414 * In particular, if there are any non-zero coefficients among
1415 * the non-parameter variables, then we take the last of these
1416 * variables. Eliminating this variable in terms of the other
1417 * variables and/or parameters does not influence the property
1418 * that all column in the initial tableau are lexicographically
1419 * positive. The row corresponding to the eliminated variable
1420 * will only have non-zero entries below the diagonal of the
1421 * initial tableau. That is, we transform
1423 * I I
1424 * 1 into a
1425 * I I
1427 * If there is no such non-parameter variable, then we are dealing with
1428 * pure parameter equality and we pick any parameter with coefficient 1 or -1
1429 * for elimination. This will ensure that the eliminated parameter
1430 * always has an integer value whenever all the other parameters are integral.
1431 * If there is no such parameter then we return -1.
1433 static int last_var_col_or_int_par_col(struct isl_tab *tab, int row)
1435 unsigned off = 2 + tab->M;
1436 int i;
1438 for (i = tab->n_var - tab->n_div - 1; i >= 0 && i >= tab->n_param; --i) {
1439 int col;
1440 if (tab->var[i].is_row)
1441 continue;
1442 col = tab->var[i].index;
1443 if (col <= tab->n_dead)
1444 continue;
1445 if (!isl_int_is_zero(tab->mat->row[row][off + col]))
1446 return col;
1448 for (i = tab->n_dead; i < tab->n_col; ++i) {
1449 if (isl_int_is_one(tab->mat->row[row][off + i]))
1450 return i;
1451 if (isl_int_is_negone(tab->mat->row[row][off + i]))
1452 return i;
1454 return -1;
1457 /* Add an equality that is known to be valid to the tableau.
1458 * We first check if we can eliminate a variable or a parameter.
1459 * If not, we add the equality as two inequalities.
1460 * In this case, the equality was a pure parameter equality and there
1461 * is no need to resolve any constraint violations.
1463 * This function assumes that at least two more rows and at least
1464 * two more elements in the constraint array are available in the tableau.
1466 static struct isl_tab *add_lexmin_valid_eq(struct isl_tab *tab, isl_int *eq)
1468 int i;
1469 int r;
1471 if (!tab)
1472 return NULL;
1473 r = isl_tab_add_row(tab, eq);
1474 if (r < 0)
1475 goto error;
1477 r = tab->con[r].index;
1478 i = last_var_col_or_int_par_col(tab, r);
1479 if (i < 0) {
1480 tab->con[r].is_nonneg = 1;
1481 if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r]) < 0)
1482 goto error;
1483 isl_seq_neg(eq, eq, 1 + tab->n_var);
1484 r = isl_tab_add_row(tab, eq);
1485 if (r < 0)
1486 goto error;
1487 tab->con[r].is_nonneg = 1;
1488 if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r]) < 0)
1489 goto error;
1490 } else {
1491 if (isl_tab_pivot(tab, r, i) < 0)
1492 goto error;
1493 if (isl_tab_kill_col(tab, i) < 0)
1494 goto error;
1495 tab->n_eq++;
1498 return tab;
1499 error:
1500 isl_tab_free(tab);
1501 return NULL;
1504 /* Check if the given row is a pure constant.
1506 static int is_constant(struct isl_tab *tab, int row)
1508 unsigned off = 2 + tab->M;
1510 return isl_seq_first_non_zero(tab->mat->row[row] + off + tab->n_dead,
1511 tab->n_col - tab->n_dead) == -1;
1514 /* Add an equality that may or may not be valid to the tableau.
1515 * If the resulting row is a pure constant, then it must be zero.
1516 * Otherwise, the resulting tableau is empty.
1518 * If the row is not a pure constant, then we add two inequalities,
1519 * each time checking that they can be satisfied.
1520 * In the end we try to use one of the two constraints to eliminate
1521 * a column.
1523 * This function assumes that at least two more rows and at least
1524 * two more elements in the constraint array are available in the tableau.
1526 static int add_lexmin_eq(struct isl_tab *tab, isl_int *eq) WARN_UNUSED;
1527 static int add_lexmin_eq(struct isl_tab *tab, isl_int *eq)
1529 int r1, r2;
1530 int row;
1531 struct isl_tab_undo *snap;
1533 if (!tab)
1534 return -1;
1535 snap = isl_tab_snap(tab);
1536 r1 = isl_tab_add_row(tab, eq);
1537 if (r1 < 0)
1538 return -1;
1539 tab->con[r1].is_nonneg = 1;
1540 if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r1]) < 0)
1541 return -1;
1543 row = tab->con[r1].index;
1544 if (is_constant(tab, row)) {
1545 if (!isl_int_is_zero(tab->mat->row[row][1]) ||
1546 (tab->M && !isl_int_is_zero(tab->mat->row[row][2]))) {
1547 if (isl_tab_mark_empty(tab) < 0)
1548 return -1;
1549 return 0;
1551 if (isl_tab_rollback(tab, snap) < 0)
1552 return -1;
1553 return 0;
1556 if (restore_lexmin(tab) < 0)
1557 return -1;
1558 if (tab->empty)
1559 return 0;
1561 isl_seq_neg(eq, eq, 1 + tab->n_var);
1563 r2 = isl_tab_add_row(tab, eq);
1564 if (r2 < 0)
1565 return -1;
1566 tab->con[r2].is_nonneg = 1;
1567 if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r2]) < 0)
1568 return -1;
1570 if (restore_lexmin(tab) < 0)
1571 return -1;
1572 if (tab->empty)
1573 return 0;
1575 if (!tab->con[r1].is_row) {
1576 if (isl_tab_kill_col(tab, tab->con[r1].index) < 0)
1577 return -1;
1578 } else if (!tab->con[r2].is_row) {
1579 if (isl_tab_kill_col(tab, tab->con[r2].index) < 0)
1580 return -1;
1583 if (tab->bmap) {
1584 tab->bmap = isl_basic_map_add_ineq(tab->bmap, eq);
1585 if (isl_tab_push(tab, isl_tab_undo_bmap_ineq) < 0)
1586 return -1;
1587 isl_seq_neg(eq, eq, 1 + tab->n_var);
1588 tab->bmap = isl_basic_map_add_ineq(tab->bmap, eq);
1589 isl_seq_neg(eq, eq, 1 + tab->n_var);
1590 if (isl_tab_push(tab, isl_tab_undo_bmap_ineq) < 0)
1591 return -1;
1592 if (!tab->bmap)
1593 return -1;
1596 return 0;
1599 /* Add an inequality to the tableau, resolving violations using
1600 * restore_lexmin.
1602 * This function assumes that at least one more row and at least
1603 * one more element in the constraint array are available in the tableau.
1605 static struct isl_tab *add_lexmin_ineq(struct isl_tab *tab, isl_int *ineq)
1607 int r;
1609 if (!tab)
1610 return NULL;
1611 if (tab->bmap) {
1612 tab->bmap = isl_basic_map_add_ineq(tab->bmap, ineq);
1613 if (isl_tab_push(tab, isl_tab_undo_bmap_ineq) < 0)
1614 goto error;
1615 if (!tab->bmap)
1616 goto error;
1618 r = isl_tab_add_row(tab, ineq);
1619 if (r < 0)
1620 goto error;
1621 tab->con[r].is_nonneg = 1;
1622 if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r]) < 0)
1623 goto error;
1624 if (isl_tab_row_is_redundant(tab, tab->con[r].index)) {
1625 if (isl_tab_mark_redundant(tab, tab->con[r].index) < 0)
1626 goto error;
1627 return tab;
1630 if (restore_lexmin(tab) < 0)
1631 goto error;
1632 if (!tab->empty && tab->con[r].is_row &&
1633 isl_tab_row_is_redundant(tab, tab->con[r].index))
1634 if (isl_tab_mark_redundant(tab, tab->con[r].index) < 0)
1635 goto error;
1636 return tab;
1637 error:
1638 isl_tab_free(tab);
1639 return NULL;
1642 /* Check if the coefficients of the parameters are all integral.
1644 static int integer_parameter(struct isl_tab *tab, int row)
1646 int i;
1647 int col;
1648 unsigned off = 2 + tab->M;
1650 for (i = 0; i < tab->n_param; ++i) {
1651 /* Eliminated parameter */
1652 if (tab->var[i].is_row)
1653 continue;
1654 col = tab->var[i].index;
1655 if (!isl_int_is_divisible_by(tab->mat->row[row][off + col],
1656 tab->mat->row[row][0]))
1657 return 0;
1659 for (i = 0; i < tab->n_div; ++i) {
1660 if (tab->var[tab->n_var - tab->n_div + i].is_row)
1661 continue;
1662 col = tab->var[tab->n_var - tab->n_div + i].index;
1663 if (!isl_int_is_divisible_by(tab->mat->row[row][off + col],
1664 tab->mat->row[row][0]))
1665 return 0;
1667 return 1;
1670 /* Check if the coefficients of the non-parameter variables are all integral.
1672 static int integer_variable(struct isl_tab *tab, int row)
1674 int i;
1675 unsigned off = 2 + tab->M;
1677 for (i = tab->n_dead; i < tab->n_col; ++i) {
1678 if (tab->col_var[i] >= 0 &&
1679 (tab->col_var[i] < tab->n_param ||
1680 tab->col_var[i] >= tab->n_var - tab->n_div))
1681 continue;
1682 if (!isl_int_is_divisible_by(tab->mat->row[row][off + i],
1683 tab->mat->row[row][0]))
1684 return 0;
1686 return 1;
1689 /* Check if the constant term is integral.
1691 static int integer_constant(struct isl_tab *tab, int row)
1693 return isl_int_is_divisible_by(tab->mat->row[row][1],
1694 tab->mat->row[row][0]);
1697 #define I_CST 1 << 0
1698 #define I_PAR 1 << 1
1699 #define I_VAR 1 << 2
1701 /* Check for next (non-parameter) variable after "var" (first if var == -1)
1702 * that is non-integer and therefore requires a cut and return
1703 * the index of the variable.
1704 * For parametric tableaus, there are three parts in a row,
1705 * the constant, the coefficients of the parameters and the rest.
1706 * For each part, we check whether the coefficients in that part
1707 * are all integral and if so, set the corresponding flag in *f.
1708 * If the constant and the parameter part are integral, then the
1709 * current sample value is integral and no cut is required
1710 * (irrespective of whether the variable part is integral).
1712 static int next_non_integer_var(struct isl_tab *tab, int var, int *f)
1714 var = var < 0 ? tab->n_param : var + 1;
1716 for (; var < tab->n_var - tab->n_div; ++var) {
1717 int flags = 0;
1718 int row;
1719 if (!tab->var[var].is_row)
1720 continue;
1721 row = tab->var[var].index;
1722 if (integer_constant(tab, row))
1723 ISL_FL_SET(flags, I_CST);
1724 if (integer_parameter(tab, row))
1725 ISL_FL_SET(flags, I_PAR);
1726 if (ISL_FL_ISSET(flags, I_CST) && ISL_FL_ISSET(flags, I_PAR))
1727 continue;
1728 if (integer_variable(tab, row))
1729 ISL_FL_SET(flags, I_VAR);
1730 *f = flags;
1731 return var;
1733 return -1;
1736 /* Check for first (non-parameter) variable that is non-integer and
1737 * therefore requires a cut and return the corresponding row.
1738 * For parametric tableaus, there are three parts in a row,
1739 * the constant, the coefficients of the parameters and the rest.
1740 * For each part, we check whether the coefficients in that part
1741 * are all integral and if so, set the corresponding flag in *f.
1742 * If the constant and the parameter part are integral, then the
1743 * current sample value is integral and no cut is required
1744 * (irrespective of whether the variable part is integral).
1746 static int first_non_integer_row(struct isl_tab *tab, int *f)
1748 int var = next_non_integer_var(tab, -1, f);
1750 return var < 0 ? -1 : tab->var[var].index;
1753 /* Add a (non-parametric) cut to cut away the non-integral sample
1754 * value of the given row.
1756 * If the row is given by
1758 * m r = f + \sum_i a_i y_i
1760 * then the cut is
1762 * c = - {-f/m} + \sum_i {a_i/m} y_i >= 0
1764 * The big parameter, if any, is ignored, since it is assumed to be big
1765 * enough to be divisible by any integer.
1766 * If the tableau is actually a parametric tableau, then this function
1767 * is only called when all coefficients of the parameters are integral.
1768 * The cut therefore has zero coefficients for the parameters.
1770 * The current value is known to be negative, so row_sign, if it
1771 * exists, is set accordingly.
1773 * Return the row of the cut or -1.
1775 static int add_cut(struct isl_tab *tab, int row)
1777 int i;
1778 int r;
1779 isl_int *r_row;
1780 unsigned off = 2 + tab->M;
1782 if (isl_tab_extend_cons(tab, 1) < 0)
1783 return -1;
1784 r = isl_tab_allocate_con(tab);
1785 if (r < 0)
1786 return -1;
1788 r_row = tab->mat->row[tab->con[r].index];
1789 isl_int_set(r_row[0], tab->mat->row[row][0]);
1790 isl_int_neg(r_row[1], tab->mat->row[row][1]);
1791 isl_int_fdiv_r(r_row[1], r_row[1], tab->mat->row[row][0]);
1792 isl_int_neg(r_row[1], r_row[1]);
1793 if (tab->M)
1794 isl_int_set_si(r_row[2], 0);
1795 for (i = 0; i < tab->n_col; ++i)
1796 isl_int_fdiv_r(r_row[off + i],
1797 tab->mat->row[row][off + i], tab->mat->row[row][0]);
1799 tab->con[r].is_nonneg = 1;
1800 if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r]) < 0)
1801 return -1;
1802 if (tab->row_sign)
1803 tab->row_sign[tab->con[r].index] = isl_tab_row_neg;
1805 return tab->con[r].index;
1808 #define CUT_ALL 1
1809 #define CUT_ONE 0
1811 /* Given a non-parametric tableau, add cuts until an integer
1812 * sample point is obtained or until the tableau is determined
1813 * to be integer infeasible.
1814 * As long as there is any non-integer value in the sample point,
1815 * we add appropriate cuts, if possible, for each of these
1816 * non-integer values and then resolve the violated
1817 * cut constraints using restore_lexmin.
1818 * If one of the corresponding rows is equal to an integral
1819 * combination of variables/constraints plus a non-integral constant,
1820 * then there is no way to obtain an integer point and we return
1821 * a tableau that is marked empty.
1822 * The parameter cutting_strategy controls the strategy used when adding cuts
1823 * to remove non-integer points. CUT_ALL adds all possible cuts
1824 * before continuing the search. CUT_ONE adds only one cut at a time.
1826 static struct isl_tab *cut_to_integer_lexmin(struct isl_tab *tab,
1827 int cutting_strategy)
1829 int var;
1830 int row;
1831 int flags;
1833 if (!tab)
1834 return NULL;
1835 if (tab->empty)
1836 return tab;
1838 while ((var = next_non_integer_var(tab, -1, &flags)) != -1) {
1839 do {
1840 if (ISL_FL_ISSET(flags, I_VAR)) {
1841 if (isl_tab_mark_empty(tab) < 0)
1842 goto error;
1843 return tab;
1845 row = tab->var[var].index;
1846 row = add_cut(tab, row);
1847 if (row < 0)
1848 goto error;
1849 if (cutting_strategy == CUT_ONE)
1850 break;
1851 } while ((var = next_non_integer_var(tab, var, &flags)) != -1);
1852 if (restore_lexmin(tab) < 0)
1853 goto error;
1854 if (tab->empty)
1855 break;
1857 return tab;
1858 error:
1859 isl_tab_free(tab);
1860 return NULL;
1863 /* Check whether all the currently active samples also satisfy the inequality
1864 * "ineq" (treated as an equality if eq is set).
1865 * Remove those samples that do not.
1867 static struct isl_tab *check_samples(struct isl_tab *tab, isl_int *ineq, int eq)
1869 int i;
1870 isl_int v;
1872 if (!tab)
1873 return NULL;
1875 isl_assert(tab->mat->ctx, tab->bmap, goto error);
1876 isl_assert(tab->mat->ctx, tab->samples, goto error);
1877 isl_assert(tab->mat->ctx, tab->samples->n_col == 1 + tab->n_var, goto error);
1879 isl_int_init(v);
1880 for (i = tab->n_outside; i < tab->n_sample; ++i) {
1881 int sgn;
1882 isl_seq_inner_product(ineq, tab->samples->row[i],
1883 1 + tab->n_var, &v);
1884 sgn = isl_int_sgn(v);
1885 if (eq ? (sgn == 0) : (sgn >= 0))
1886 continue;
1887 tab = isl_tab_drop_sample(tab, i);
1888 if (!tab)
1889 break;
1891 isl_int_clear(v);
1893 return tab;
1894 error:
1895 isl_tab_free(tab);
1896 return NULL;
1899 /* Check whether the sample value of the tableau is finite,
1900 * i.e., either the tableau does not use a big parameter, or
1901 * all values of the variables are equal to the big parameter plus
1902 * some constant. This constant is the actual sample value.
1904 static int sample_is_finite(struct isl_tab *tab)
1906 int i;
1908 if (!tab->M)
1909 return 1;
1911 for (i = 0; i < tab->n_var; ++i) {
1912 int row;
1913 if (!tab->var[i].is_row)
1914 return 0;
1915 row = tab->var[i].index;
1916 if (isl_int_ne(tab->mat->row[row][0], tab->mat->row[row][2]))
1917 return 0;
1919 return 1;
1922 /* Check if the context tableau of sol has any integer points.
1923 * Leave tab in empty state if no integer point can be found.
1924 * If an integer point can be found and if moreover it is finite,
1925 * then it is added to the list of sample values.
1927 * This function is only called when none of the currently active sample
1928 * values satisfies the most recently added constraint.
1930 static struct isl_tab *check_integer_feasible(struct isl_tab *tab)
1932 struct isl_tab_undo *snap;
1934 if (!tab)
1935 return NULL;
1937 snap = isl_tab_snap(tab);
1938 if (isl_tab_push_basis(tab) < 0)
1939 goto error;
1941 tab = cut_to_integer_lexmin(tab, CUT_ALL);
1942 if (!tab)
1943 goto error;
1945 if (!tab->empty && sample_is_finite(tab)) {
1946 struct isl_vec *sample;
1948 sample = isl_tab_get_sample_value(tab);
1950 if (isl_tab_add_sample(tab, sample) < 0)
1951 goto error;
1954 if (!tab->empty && isl_tab_rollback(tab, snap) < 0)
1955 goto error;
1957 return tab;
1958 error:
1959 isl_tab_free(tab);
1960 return NULL;
1963 /* Check if any of the currently active sample values satisfies
1964 * the inequality "ineq" (an equality if eq is set).
1966 static int tab_has_valid_sample(struct isl_tab *tab, isl_int *ineq, int eq)
1968 int i;
1969 isl_int v;
1971 if (!tab)
1972 return -1;
1974 isl_assert(tab->mat->ctx, tab->bmap, return -1);
1975 isl_assert(tab->mat->ctx, tab->samples, return -1);
1976 isl_assert(tab->mat->ctx, tab->samples->n_col == 1 + tab->n_var, return -1);
1978 isl_int_init(v);
1979 for (i = tab->n_outside; i < tab->n_sample; ++i) {
1980 int sgn;
1981 isl_seq_inner_product(ineq, tab->samples->row[i],
1982 1 + tab->n_var, &v);
1983 sgn = isl_int_sgn(v);
1984 if (eq ? (sgn == 0) : (sgn >= 0))
1985 break;
1987 isl_int_clear(v);
1989 return i < tab->n_sample;
1992 /* Insert a div specified by "div" to the tableau "tab" at position "pos" and
1993 * return isl_bool_true if the div is obviously non-negative.
1995 static isl_bool context_tab_insert_div(struct isl_tab *tab, int pos,
1996 __isl_keep isl_vec *div,
1997 isl_stat (*add_ineq)(void *user, isl_int *), void *user)
1999 int i;
2000 int r;
2001 struct isl_mat *samples;
2002 int nonneg;
2004 r = isl_tab_insert_div(tab, pos, div, add_ineq, user);
2005 if (r < 0)
2006 return isl_bool_error;
2007 nonneg = tab->var[r].is_nonneg;
2008 tab->var[r].frozen = 1;
2010 samples = isl_mat_extend(tab->samples,
2011 tab->n_sample, 1 + tab->n_var);
2012 tab->samples = samples;
2013 if (!samples)
2014 return isl_bool_error;
2015 for (i = tab->n_outside; i < samples->n_row; ++i) {
2016 isl_seq_inner_product(div->el + 1, samples->row[i],
2017 div->size - 1, &samples->row[i][samples->n_col - 1]);
2018 isl_int_fdiv_q(samples->row[i][samples->n_col - 1],
2019 samples->row[i][samples->n_col - 1], div->el[0]);
2021 tab->samples = isl_mat_move_cols(tab->samples, 1 + pos,
2022 1 + tab->n_var - 1, 1);
2023 if (!tab->samples)
2024 return isl_bool_error;
2026 return nonneg;
2029 /* Add a div specified by "div" to both the main tableau and
2030 * the context tableau. In case of the main tableau, we only
2031 * need to add an extra div. In the context tableau, we also
2032 * need to express the meaning of the div.
2033 * Return the index of the div or -1 if anything went wrong.
2035 * The new integer division is added before any unknown integer
2036 * divisions in the context to ensure that it does not get
2037 * equated to some linear combination involving unknown integer
2038 * divisions.
2040 static int add_div(struct isl_tab *tab, struct isl_context *context,
2041 __isl_keep isl_vec *div)
2043 int r;
2044 int pos;
2045 isl_bool nonneg;
2046 struct isl_tab *context_tab = context->op->peek_tab(context);
2048 if (!tab || !context_tab)
2049 goto error;
2051 pos = context_tab->n_var - context->n_unknown;
2052 if ((nonneg = context->op->insert_div(context, pos, div)) < 0)
2053 goto error;
2055 if (!context->op->is_ok(context))
2056 goto error;
2058 pos = tab->n_var - context->n_unknown;
2059 if (isl_tab_extend_vars(tab, 1) < 0)
2060 goto error;
2061 r = isl_tab_insert_var(tab, pos);
2062 if (r < 0)
2063 goto error;
2064 if (nonneg)
2065 tab->var[r].is_nonneg = 1;
2066 tab->var[r].frozen = 1;
2067 tab->n_div++;
2069 return tab->n_div - 1 - context->n_unknown;
2070 error:
2071 context->op->invalidate(context);
2072 return -1;
2075 static int find_div(struct isl_tab *tab, isl_int *div, isl_int denom)
2077 int i;
2078 unsigned total = isl_basic_map_total_dim(tab->bmap);
2080 for (i = 0; i < tab->bmap->n_div; ++i) {
2081 if (isl_int_ne(tab->bmap->div[i][0], denom))
2082 continue;
2083 if (!isl_seq_eq(tab->bmap->div[i] + 1, div, 1 + total))
2084 continue;
2085 return i;
2087 return -1;
2090 /* Return the index of a div that corresponds to "div".
2091 * We first check if we already have such a div and if not, we create one.
2093 static int get_div(struct isl_tab *tab, struct isl_context *context,
2094 struct isl_vec *div)
2096 int d;
2097 struct isl_tab *context_tab = context->op->peek_tab(context);
2099 if (!context_tab)
2100 return -1;
2102 d = find_div(context_tab, div->el + 1, div->el[0]);
2103 if (d != -1)
2104 return d;
2106 return add_div(tab, context, div);
2109 /* Add a parametric cut to cut away the non-integral sample value
2110 * of the give row.
2111 * Let a_i be the coefficients of the constant term and the parameters
2112 * and let b_i be the coefficients of the variables or constraints
2113 * in basis of the tableau.
2114 * Let q be the div q = floor(\sum_i {-a_i} y_i).
2116 * The cut is expressed as
2118 * c = \sum_i -{-a_i} y_i + \sum_i {b_i} x_i + q >= 0
2120 * If q did not already exist in the context tableau, then it is added first.
2121 * If q is in a column of the main tableau then the "+ q" can be accomplished
2122 * by setting the corresponding entry to the denominator of the constraint.
2123 * If q happens to be in a row of the main tableau, then the corresponding
2124 * row needs to be added instead (taking care of the denominators).
2125 * Note that this is very unlikely, but perhaps not entirely impossible.
2127 * The current value of the cut is known to be negative (or at least
2128 * non-positive), so row_sign is set accordingly.
2130 * Return the row of the cut or -1.
2132 static int add_parametric_cut(struct isl_tab *tab, int row,
2133 struct isl_context *context)
2135 struct isl_vec *div;
2136 int d;
2137 int i;
2138 int r;
2139 isl_int *r_row;
2140 int col;
2141 int n;
2142 unsigned off = 2 + tab->M;
2144 if (!context)
2145 return -1;
2147 div = get_row_parameter_div(tab, row);
2148 if (!div)
2149 return -1;
2151 n = tab->n_div - context->n_unknown;
2152 d = context->op->get_div(context, tab, div);
2153 isl_vec_free(div);
2154 if (d < 0)
2155 return -1;
2157 if (isl_tab_extend_cons(tab, 1) < 0)
2158 return -1;
2159 r = isl_tab_allocate_con(tab);
2160 if (r < 0)
2161 return -1;
2163 r_row = tab->mat->row[tab->con[r].index];
2164 isl_int_set(r_row[0], tab->mat->row[row][0]);
2165 isl_int_neg(r_row[1], tab->mat->row[row][1]);
2166 isl_int_fdiv_r(r_row[1], r_row[1], tab->mat->row[row][0]);
2167 isl_int_neg(r_row[1], r_row[1]);
2168 if (tab->M)
2169 isl_int_set_si(r_row[2], 0);
2170 for (i = 0; i < tab->n_param; ++i) {
2171 if (tab->var[i].is_row)
2172 continue;
2173 col = tab->var[i].index;
2174 isl_int_neg(r_row[off + col], tab->mat->row[row][off + col]);
2175 isl_int_fdiv_r(r_row[off + col], r_row[off + col],
2176 tab->mat->row[row][0]);
2177 isl_int_neg(r_row[off + col], r_row[off + col]);
2179 for (i = 0; i < tab->n_div; ++i) {
2180 if (tab->var[tab->n_var - tab->n_div + i].is_row)
2181 continue;
2182 col = tab->var[tab->n_var - tab->n_div + i].index;
2183 isl_int_neg(r_row[off + col], tab->mat->row[row][off + col]);
2184 isl_int_fdiv_r(r_row[off + col], r_row[off + col],
2185 tab->mat->row[row][0]);
2186 isl_int_neg(r_row[off + col], r_row[off + col]);
2188 for (i = 0; i < tab->n_col; ++i) {
2189 if (tab->col_var[i] >= 0 &&
2190 (tab->col_var[i] < tab->n_param ||
2191 tab->col_var[i] >= tab->n_var - tab->n_div))
2192 continue;
2193 isl_int_fdiv_r(r_row[off + i],
2194 tab->mat->row[row][off + i], tab->mat->row[row][0]);
2196 if (tab->var[tab->n_var - tab->n_div + d].is_row) {
2197 isl_int gcd;
2198 int d_row = tab->var[tab->n_var - tab->n_div + d].index;
2199 isl_int_init(gcd);
2200 isl_int_gcd(gcd, tab->mat->row[d_row][0], r_row[0]);
2201 isl_int_divexact(r_row[0], r_row[0], gcd);
2202 isl_int_divexact(gcd, tab->mat->row[d_row][0], gcd);
2203 isl_seq_combine(r_row + 1, gcd, r_row + 1,
2204 r_row[0], tab->mat->row[d_row] + 1,
2205 off - 1 + tab->n_col);
2206 isl_int_mul(r_row[0], r_row[0], tab->mat->row[d_row][0]);
2207 isl_int_clear(gcd);
2208 } else {
2209 col = tab->var[tab->n_var - tab->n_div + d].index;
2210 isl_int_set(r_row[off + col], tab->mat->row[row][0]);
2213 tab->con[r].is_nonneg = 1;
2214 if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r]) < 0)
2215 return -1;
2216 if (tab->row_sign)
2217 tab->row_sign[tab->con[r].index] = isl_tab_row_neg;
2219 row = tab->con[r].index;
2221 if (d >= n && context->op->detect_equalities(context, tab) < 0)
2222 return -1;
2224 return row;
2227 /* Construct a tableau for bmap that can be used for computing
2228 * the lexicographic minimum (or maximum) of bmap.
2229 * If not NULL, then dom is the domain where the minimum
2230 * should be computed. In this case, we set up a parametric
2231 * tableau with row signs (initialized to "unknown").
2232 * If M is set, then the tableau will use a big parameter.
2233 * If max is set, then a maximum should be computed instead of a minimum.
2234 * This means that for each variable x, the tableau will contain the variable
2235 * x' = M - x, rather than x' = M + x. This in turn means that the coefficient
2236 * of the variables in all constraints are negated prior to adding them
2237 * to the tableau.
2239 static __isl_give struct isl_tab *tab_for_lexmin(__isl_keep isl_basic_map *bmap,
2240 __isl_keep isl_basic_set *dom, unsigned M, int max)
2242 int i;
2243 struct isl_tab *tab;
2244 unsigned n_var;
2245 unsigned o_var;
2247 tab = isl_tab_alloc(bmap->ctx, 2 * bmap->n_eq + bmap->n_ineq + 1,
2248 isl_basic_map_total_dim(bmap), M);
2249 if (!tab)
2250 return NULL;
2252 tab->rational = ISL_F_ISSET(bmap, ISL_BASIC_MAP_RATIONAL);
2253 if (dom) {
2254 tab->n_param = isl_basic_set_total_dim(dom) - dom->n_div;
2255 tab->n_div = dom->n_div;
2256 tab->row_sign = isl_calloc_array(bmap->ctx,
2257 enum isl_tab_row_sign, tab->mat->n_row);
2258 if (tab->mat->n_row && !tab->row_sign)
2259 goto error;
2261 if (ISL_F_ISSET(bmap, ISL_BASIC_MAP_EMPTY)) {
2262 if (isl_tab_mark_empty(tab) < 0)
2263 goto error;
2264 return tab;
2267 for (i = tab->n_param; i < tab->n_var - tab->n_div; ++i) {
2268 tab->var[i].is_nonneg = 1;
2269 tab->var[i].frozen = 1;
2271 o_var = 1 + tab->n_param;
2272 n_var = tab->n_var - tab->n_param - tab->n_div;
2273 for (i = 0; i < bmap->n_eq; ++i) {
2274 if (max)
2275 isl_seq_neg(bmap->eq[i] + o_var,
2276 bmap->eq[i] + o_var, n_var);
2277 tab = add_lexmin_valid_eq(tab, bmap->eq[i]);
2278 if (max)
2279 isl_seq_neg(bmap->eq[i] + o_var,
2280 bmap->eq[i] + o_var, n_var);
2281 if (!tab || tab->empty)
2282 return tab;
2284 if (bmap->n_eq && restore_lexmin(tab) < 0)
2285 goto error;
2286 for (i = 0; i < bmap->n_ineq; ++i) {
2287 if (max)
2288 isl_seq_neg(bmap->ineq[i] + o_var,
2289 bmap->ineq[i] + o_var, n_var);
2290 tab = add_lexmin_ineq(tab, bmap->ineq[i]);
2291 if (max)
2292 isl_seq_neg(bmap->ineq[i] + o_var,
2293 bmap->ineq[i] + o_var, n_var);
2294 if (!tab || tab->empty)
2295 return tab;
2297 return tab;
2298 error:
2299 isl_tab_free(tab);
2300 return NULL;
2303 /* Given a main tableau where more than one row requires a split,
2304 * determine and return the "best" row to split on.
2306 * Given two rows in the main tableau, if the inequality corresponding
2307 * to the first row is redundant with respect to that of the second row
2308 * in the current tableau, then it is better to split on the second row,
2309 * since in the positive part, both rows will be positive.
2310 * (In the negative part a pivot will have to be performed and just about
2311 * anything can happen to the sign of the other row.)
2313 * As a simple heuristic, we therefore select the row that makes the most
2314 * of the other rows redundant.
2316 * Perhaps it would also be useful to look at the number of constraints
2317 * that conflict with any given constraint.
2319 * best is the best row so far (-1 when we have not found any row yet).
2320 * best_r is the number of other rows made redundant by row best.
2321 * When best is still -1, bset_r is meaningless, but it is initialized
2322 * to some arbitrary value (0) anyway. Without this redundant initialization
2323 * valgrind may warn about uninitialized memory accesses when isl
2324 * is compiled with some versions of gcc.
2326 static int best_split(struct isl_tab *tab, struct isl_tab *context_tab)
2328 struct isl_tab_undo *snap;
2329 int split;
2330 int row;
2331 int best = -1;
2332 int best_r = 0;
2334 if (isl_tab_extend_cons(context_tab, 2) < 0)
2335 return -1;
2337 snap = isl_tab_snap(context_tab);
2339 for (split = tab->n_redundant; split < tab->n_row; ++split) {
2340 struct isl_tab_undo *snap2;
2341 struct isl_vec *ineq = NULL;
2342 int r = 0;
2343 int ok;
2345 if (!isl_tab_var_from_row(tab, split)->is_nonneg)
2346 continue;
2347 if (tab->row_sign[split] != isl_tab_row_any)
2348 continue;
2350 ineq = get_row_parameter_ineq(tab, split);
2351 if (!ineq)
2352 return -1;
2353 ok = isl_tab_add_ineq(context_tab, ineq->el) >= 0;
2354 isl_vec_free(ineq);
2355 if (!ok)
2356 return -1;
2358 snap2 = isl_tab_snap(context_tab);
2360 for (row = tab->n_redundant; row < tab->n_row; ++row) {
2361 struct isl_tab_var *var;
2363 if (row == split)
2364 continue;
2365 if (!isl_tab_var_from_row(tab, row)->is_nonneg)
2366 continue;
2367 if (tab->row_sign[row] != isl_tab_row_any)
2368 continue;
2370 ineq = get_row_parameter_ineq(tab, row);
2371 if (!ineq)
2372 return -1;
2373 ok = isl_tab_add_ineq(context_tab, ineq->el) >= 0;
2374 isl_vec_free(ineq);
2375 if (!ok)
2376 return -1;
2377 var = &context_tab->con[context_tab->n_con - 1];
2378 if (!context_tab->empty &&
2379 !isl_tab_min_at_most_neg_one(context_tab, var))
2380 r++;
2381 if (isl_tab_rollback(context_tab, snap2) < 0)
2382 return -1;
2384 if (best == -1 || r > best_r) {
2385 best = split;
2386 best_r = r;
2388 if (isl_tab_rollback(context_tab, snap) < 0)
2389 return -1;
2392 return best;
2395 static struct isl_basic_set *context_lex_peek_basic_set(
2396 struct isl_context *context)
2398 struct isl_context_lex *clex = (struct isl_context_lex *)context;
2399 if (!clex->tab)
2400 return NULL;
2401 return isl_tab_peek_bset(clex->tab);
2404 static struct isl_tab *context_lex_peek_tab(struct isl_context *context)
2406 struct isl_context_lex *clex = (struct isl_context_lex *)context;
2407 return clex->tab;
2410 static void context_lex_add_eq(struct isl_context *context, isl_int *eq,
2411 int check, int update)
2413 struct isl_context_lex *clex = (struct isl_context_lex *)context;
2414 if (isl_tab_extend_cons(clex->tab, 2) < 0)
2415 goto error;
2416 if (add_lexmin_eq(clex->tab, eq) < 0)
2417 goto error;
2418 if (check) {
2419 int v = tab_has_valid_sample(clex->tab, eq, 1);
2420 if (v < 0)
2421 goto error;
2422 if (!v)
2423 clex->tab = check_integer_feasible(clex->tab);
2425 if (update)
2426 clex->tab = check_samples(clex->tab, eq, 1);
2427 return;
2428 error:
2429 isl_tab_free(clex->tab);
2430 clex->tab = NULL;
2433 static void context_lex_add_ineq(struct isl_context *context, isl_int *ineq,
2434 int check, int update)
2436 struct isl_context_lex *clex = (struct isl_context_lex *)context;
2437 if (isl_tab_extend_cons(clex->tab, 1) < 0)
2438 goto error;
2439 clex->tab = add_lexmin_ineq(clex->tab, ineq);
2440 if (check) {
2441 int v = tab_has_valid_sample(clex->tab, ineq, 0);
2442 if (v < 0)
2443 goto error;
2444 if (!v)
2445 clex->tab = check_integer_feasible(clex->tab);
2447 if (update)
2448 clex->tab = check_samples(clex->tab, ineq, 0);
2449 return;
2450 error:
2451 isl_tab_free(clex->tab);
2452 clex->tab = NULL;
2455 static isl_stat context_lex_add_ineq_wrap(void *user, isl_int *ineq)
2457 struct isl_context *context = (struct isl_context *)user;
2458 context_lex_add_ineq(context, ineq, 0, 0);
2459 return context->op->is_ok(context) ? isl_stat_ok : isl_stat_error;
2462 /* Check which signs can be obtained by "ineq" on all the currently
2463 * active sample values. See row_sign for more information.
2465 static enum isl_tab_row_sign tab_ineq_sign(struct isl_tab *tab, isl_int *ineq,
2466 int strict)
2468 int i;
2469 int sgn;
2470 isl_int tmp;
2471 enum isl_tab_row_sign res = isl_tab_row_unknown;
2473 isl_assert(tab->mat->ctx, tab->samples, return isl_tab_row_unknown);
2474 isl_assert(tab->mat->ctx, tab->samples->n_col == 1 + tab->n_var,
2475 return isl_tab_row_unknown);
2477 isl_int_init(tmp);
2478 for (i = tab->n_outside; i < tab->n_sample; ++i) {
2479 isl_seq_inner_product(tab->samples->row[i], ineq,
2480 1 + tab->n_var, &tmp);
2481 sgn = isl_int_sgn(tmp);
2482 if (sgn > 0 || (sgn == 0 && strict)) {
2483 if (res == isl_tab_row_unknown)
2484 res = isl_tab_row_pos;
2485 if (res == isl_tab_row_neg)
2486 res = isl_tab_row_any;
2488 if (sgn < 0) {
2489 if (res == isl_tab_row_unknown)
2490 res = isl_tab_row_neg;
2491 if (res == isl_tab_row_pos)
2492 res = isl_tab_row_any;
2494 if (res == isl_tab_row_any)
2495 break;
2497 isl_int_clear(tmp);
2499 return res;
2502 static enum isl_tab_row_sign context_lex_ineq_sign(struct isl_context *context,
2503 isl_int *ineq, int strict)
2505 struct isl_context_lex *clex = (struct isl_context_lex *)context;
2506 return tab_ineq_sign(clex->tab, ineq, strict);
2509 /* Check whether "ineq" can be added to the tableau without rendering
2510 * it infeasible.
2512 static int context_lex_test_ineq(struct isl_context *context, isl_int *ineq)
2514 struct isl_context_lex *clex = (struct isl_context_lex *)context;
2515 struct isl_tab_undo *snap;
2516 int feasible;
2518 if (!clex->tab)
2519 return -1;
2521 if (isl_tab_extend_cons(clex->tab, 1) < 0)
2522 return -1;
2524 snap = isl_tab_snap(clex->tab);
2525 if (isl_tab_push_basis(clex->tab) < 0)
2526 return -1;
2527 clex->tab = add_lexmin_ineq(clex->tab, ineq);
2528 clex->tab = check_integer_feasible(clex->tab);
2529 if (!clex->tab)
2530 return -1;
2531 feasible = !clex->tab->empty;
2532 if (isl_tab_rollback(clex->tab, snap) < 0)
2533 return -1;
2535 return feasible;
2538 static int context_lex_get_div(struct isl_context *context, struct isl_tab *tab,
2539 struct isl_vec *div)
2541 return get_div(tab, context, div);
2544 /* Insert a div specified by "div" to the context tableau at position "pos" and
2545 * return isl_bool_true if the div is obviously non-negative.
2546 * context_tab_add_div will always return isl_bool_true, because all variables
2547 * in a isl_context_lex tableau are non-negative.
2548 * However, if we are using a big parameter in the context, then this only
2549 * reflects the non-negativity of the variable used to _encode_ the
2550 * div, i.e., div' = M + div, so we can't draw any conclusions.
2552 static isl_bool context_lex_insert_div(struct isl_context *context, int pos,
2553 __isl_keep isl_vec *div)
2555 struct isl_context_lex *clex = (struct isl_context_lex *)context;
2556 isl_bool nonneg;
2557 nonneg = context_tab_insert_div(clex->tab, pos, div,
2558 context_lex_add_ineq_wrap, context);
2559 if (nonneg < 0)
2560 return isl_bool_error;
2561 if (clex->tab->M)
2562 return isl_bool_false;
2563 return nonneg;
2566 static int context_lex_detect_equalities(struct isl_context *context,
2567 struct isl_tab *tab)
2569 return 0;
2572 static int context_lex_best_split(struct isl_context *context,
2573 struct isl_tab *tab)
2575 struct isl_context_lex *clex = (struct isl_context_lex *)context;
2576 struct isl_tab_undo *snap;
2577 int r;
2579 snap = isl_tab_snap(clex->tab);
2580 if (isl_tab_push_basis(clex->tab) < 0)
2581 return -1;
2582 r = best_split(tab, clex->tab);
2584 if (r >= 0 && isl_tab_rollback(clex->tab, snap) < 0)
2585 return -1;
2587 return r;
2590 static int context_lex_is_empty(struct isl_context *context)
2592 struct isl_context_lex *clex = (struct isl_context_lex *)context;
2593 if (!clex->tab)
2594 return -1;
2595 return clex->tab->empty;
2598 static void *context_lex_save(struct isl_context *context)
2600 struct isl_context_lex *clex = (struct isl_context_lex *)context;
2601 struct isl_tab_undo *snap;
2603 snap = isl_tab_snap(clex->tab);
2604 if (isl_tab_push_basis(clex->tab) < 0)
2605 return NULL;
2606 if (isl_tab_save_samples(clex->tab) < 0)
2607 return NULL;
2609 return snap;
2612 static void context_lex_restore(struct isl_context *context, void *save)
2614 struct isl_context_lex *clex = (struct isl_context_lex *)context;
2615 if (isl_tab_rollback(clex->tab, (struct isl_tab_undo *)save) < 0) {
2616 isl_tab_free(clex->tab);
2617 clex->tab = NULL;
2621 static void context_lex_discard(void *save)
2625 static int context_lex_is_ok(struct isl_context *context)
2627 struct isl_context_lex *clex = (struct isl_context_lex *)context;
2628 return !!clex->tab;
2631 /* For each variable in the context tableau, check if the variable can
2632 * only attain non-negative values. If so, mark the parameter as non-negative
2633 * in the main tableau. This allows for a more direct identification of some
2634 * cases of violated constraints.
2636 static struct isl_tab *tab_detect_nonnegative_parameters(struct isl_tab *tab,
2637 struct isl_tab *context_tab)
2639 int i;
2640 struct isl_tab_undo *snap;
2641 struct isl_vec *ineq = NULL;
2642 struct isl_tab_var *var;
2643 int n;
2645 if (context_tab->n_var == 0)
2646 return tab;
2648 ineq = isl_vec_alloc(tab->mat->ctx, 1 + context_tab->n_var);
2649 if (!ineq)
2650 goto error;
2652 if (isl_tab_extend_cons(context_tab, 1) < 0)
2653 goto error;
2655 snap = isl_tab_snap(context_tab);
2657 n = 0;
2658 isl_seq_clr(ineq->el, ineq->size);
2659 for (i = 0; i < context_tab->n_var; ++i) {
2660 isl_int_set_si(ineq->el[1 + i], 1);
2661 if (isl_tab_add_ineq(context_tab, ineq->el) < 0)
2662 goto error;
2663 var = &context_tab->con[context_tab->n_con - 1];
2664 if (!context_tab->empty &&
2665 !isl_tab_min_at_most_neg_one(context_tab, var)) {
2666 int j = i;
2667 if (i >= tab->n_param)
2668 j = i - tab->n_param + tab->n_var - tab->n_div;
2669 tab->var[j].is_nonneg = 1;
2670 n++;
2672 isl_int_set_si(ineq->el[1 + i], 0);
2673 if (isl_tab_rollback(context_tab, snap) < 0)
2674 goto error;
2677 if (context_tab->M && n == context_tab->n_var) {
2678 context_tab->mat = isl_mat_drop_cols(context_tab->mat, 2, 1);
2679 context_tab->M = 0;
2682 isl_vec_free(ineq);
2683 return tab;
2684 error:
2685 isl_vec_free(ineq);
2686 isl_tab_free(tab);
2687 return NULL;
2690 static struct isl_tab *context_lex_detect_nonnegative_parameters(
2691 struct isl_context *context, struct isl_tab *tab)
2693 struct isl_context_lex *clex = (struct isl_context_lex *)context;
2694 struct isl_tab_undo *snap;
2696 if (!tab)
2697 return NULL;
2699 snap = isl_tab_snap(clex->tab);
2700 if (isl_tab_push_basis(clex->tab) < 0)
2701 goto error;
2703 tab = tab_detect_nonnegative_parameters(tab, clex->tab);
2705 if (isl_tab_rollback(clex->tab, snap) < 0)
2706 goto error;
2708 return tab;
2709 error:
2710 isl_tab_free(tab);
2711 return NULL;
2714 static void context_lex_invalidate(struct isl_context *context)
2716 struct isl_context_lex *clex = (struct isl_context_lex *)context;
2717 isl_tab_free(clex->tab);
2718 clex->tab = NULL;
2721 static __isl_null struct isl_context *context_lex_free(
2722 struct isl_context *context)
2724 struct isl_context_lex *clex = (struct isl_context_lex *)context;
2725 isl_tab_free(clex->tab);
2726 free(clex);
2728 return NULL;
2731 struct isl_context_op isl_context_lex_op = {
2732 context_lex_detect_nonnegative_parameters,
2733 context_lex_peek_basic_set,
2734 context_lex_peek_tab,
2735 context_lex_add_eq,
2736 context_lex_add_ineq,
2737 context_lex_ineq_sign,
2738 context_lex_test_ineq,
2739 context_lex_get_div,
2740 context_lex_insert_div,
2741 context_lex_detect_equalities,
2742 context_lex_best_split,
2743 context_lex_is_empty,
2744 context_lex_is_ok,
2745 context_lex_save,
2746 context_lex_restore,
2747 context_lex_discard,
2748 context_lex_invalidate,
2749 context_lex_free,
2752 static struct isl_tab *context_tab_for_lexmin(__isl_take isl_basic_set *bset)
2754 struct isl_tab *tab;
2756 if (!bset)
2757 return NULL;
2758 tab = tab_for_lexmin(bset_to_bmap(bset), NULL, 1, 0);
2759 if (isl_tab_track_bset(tab, bset) < 0)
2760 goto error;
2761 tab = isl_tab_init_samples(tab);
2762 return tab;
2763 error:
2764 isl_tab_free(tab);
2765 return NULL;
2768 static struct isl_context *isl_context_lex_alloc(struct isl_basic_set *dom)
2770 struct isl_context_lex *clex;
2772 if (!dom)
2773 return NULL;
2775 clex = isl_alloc_type(dom->ctx, struct isl_context_lex);
2776 if (!clex)
2777 return NULL;
2779 clex->context.op = &isl_context_lex_op;
2781 clex->tab = context_tab_for_lexmin(isl_basic_set_copy(dom));
2782 if (restore_lexmin(clex->tab) < 0)
2783 goto error;
2784 clex->tab = check_integer_feasible(clex->tab);
2785 if (!clex->tab)
2786 goto error;
2788 return &clex->context;
2789 error:
2790 clex->context.op->free(&clex->context);
2791 return NULL;
2794 /* Representation of the context when using generalized basis reduction.
2796 * "shifted" contains the offsets of the unit hypercubes that lie inside the
2797 * context. Any rational point in "shifted" can therefore be rounded
2798 * up to an integer point in the context.
2799 * If the context is constrained by any equality, then "shifted" is not used
2800 * as it would be empty.
2802 struct isl_context_gbr {
2803 struct isl_context context;
2804 struct isl_tab *tab;
2805 struct isl_tab *shifted;
2806 struct isl_tab *cone;
2809 static struct isl_tab *context_gbr_detect_nonnegative_parameters(
2810 struct isl_context *context, struct isl_tab *tab)
2812 struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
2813 if (!tab)
2814 return NULL;
2815 return tab_detect_nonnegative_parameters(tab, cgbr->tab);
2818 static struct isl_basic_set *context_gbr_peek_basic_set(
2819 struct isl_context *context)
2821 struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
2822 if (!cgbr->tab)
2823 return NULL;
2824 return isl_tab_peek_bset(cgbr->tab);
2827 static struct isl_tab *context_gbr_peek_tab(struct isl_context *context)
2829 struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
2830 return cgbr->tab;
2833 /* Initialize the "shifted" tableau of the context, which
2834 * contains the constraints of the original tableau shifted
2835 * by the sum of all negative coefficients. This ensures
2836 * that any rational point in the shifted tableau can
2837 * be rounded up to yield an integer point in the original tableau.
2839 static void gbr_init_shifted(struct isl_context_gbr *cgbr)
2841 int i, j;
2842 struct isl_vec *cst;
2843 struct isl_basic_set *bset = isl_tab_peek_bset(cgbr->tab);
2844 unsigned dim = isl_basic_set_total_dim(bset);
2846 cst = isl_vec_alloc(cgbr->tab->mat->ctx, bset->n_ineq);
2847 if (!cst)
2848 return;
2850 for (i = 0; i < bset->n_ineq; ++i) {
2851 isl_int_set(cst->el[i], bset->ineq[i][0]);
2852 for (j = 0; j < dim; ++j) {
2853 if (!isl_int_is_neg(bset->ineq[i][1 + j]))
2854 continue;
2855 isl_int_add(bset->ineq[i][0], bset->ineq[i][0],
2856 bset->ineq[i][1 + j]);
2860 cgbr->shifted = isl_tab_from_basic_set(bset, 0);
2862 for (i = 0; i < bset->n_ineq; ++i)
2863 isl_int_set(bset->ineq[i][0], cst->el[i]);
2865 isl_vec_free(cst);
2868 /* Check if the shifted tableau is non-empty, and if so
2869 * use the sample point to construct an integer point
2870 * of the context tableau.
2872 static struct isl_vec *gbr_get_shifted_sample(struct isl_context_gbr *cgbr)
2874 struct isl_vec *sample;
2876 if (!cgbr->shifted)
2877 gbr_init_shifted(cgbr);
2878 if (!cgbr->shifted)
2879 return NULL;
2880 if (cgbr->shifted->empty)
2881 return isl_vec_alloc(cgbr->tab->mat->ctx, 0);
2883 sample = isl_tab_get_sample_value(cgbr->shifted);
2884 sample = isl_vec_ceil(sample);
2886 return sample;
2889 static __isl_give isl_basic_set *drop_constant_terms(
2890 __isl_take isl_basic_set *bset)
2892 int i;
2894 if (!bset)
2895 return NULL;
2897 for (i = 0; i < bset->n_eq; ++i)
2898 isl_int_set_si(bset->eq[i][0], 0);
2900 for (i = 0; i < bset->n_ineq; ++i)
2901 isl_int_set_si(bset->ineq[i][0], 0);
2903 return bset;
2906 static int use_shifted(struct isl_context_gbr *cgbr)
2908 if (!cgbr->tab)
2909 return 0;
2910 return cgbr->tab->bmap->n_eq == 0 && cgbr->tab->bmap->n_div == 0;
2913 static struct isl_vec *gbr_get_sample(struct isl_context_gbr *cgbr)
2915 struct isl_basic_set *bset;
2916 struct isl_basic_set *cone;
2918 if (isl_tab_sample_is_integer(cgbr->tab))
2919 return isl_tab_get_sample_value(cgbr->tab);
2921 if (use_shifted(cgbr)) {
2922 struct isl_vec *sample;
2924 sample = gbr_get_shifted_sample(cgbr);
2925 if (!sample || sample->size > 0)
2926 return sample;
2928 isl_vec_free(sample);
2931 if (!cgbr->cone) {
2932 bset = isl_tab_peek_bset(cgbr->tab);
2933 cgbr->cone = isl_tab_from_recession_cone(bset, 0);
2934 if (!cgbr->cone)
2935 return NULL;
2936 if (isl_tab_track_bset(cgbr->cone,
2937 isl_basic_set_copy(bset)) < 0)
2938 return NULL;
2940 if (isl_tab_detect_implicit_equalities(cgbr->cone) < 0)
2941 return NULL;
2943 if (cgbr->cone->n_dead == cgbr->cone->n_col) {
2944 struct isl_vec *sample;
2945 struct isl_tab_undo *snap;
2947 if (cgbr->tab->basis) {
2948 if (cgbr->tab->basis->n_col != 1 + cgbr->tab->n_var) {
2949 isl_mat_free(cgbr->tab->basis);
2950 cgbr->tab->basis = NULL;
2952 cgbr->tab->n_zero = 0;
2953 cgbr->tab->n_unbounded = 0;
2956 snap = isl_tab_snap(cgbr->tab);
2958 sample = isl_tab_sample(cgbr->tab);
2960 if (!sample || isl_tab_rollback(cgbr->tab, snap) < 0) {
2961 isl_vec_free(sample);
2962 return NULL;
2965 return sample;
2968 cone = isl_basic_set_dup(isl_tab_peek_bset(cgbr->cone));
2969 cone = drop_constant_terms(cone);
2970 cone = isl_basic_set_update_from_tab(cone, cgbr->cone);
2971 cone = isl_basic_set_underlying_set(cone);
2972 cone = isl_basic_set_gauss(cone, NULL);
2974 bset = isl_basic_set_dup(isl_tab_peek_bset(cgbr->tab));
2975 bset = isl_basic_set_update_from_tab(bset, cgbr->tab);
2976 bset = isl_basic_set_underlying_set(bset);
2977 bset = isl_basic_set_gauss(bset, NULL);
2979 return isl_basic_set_sample_with_cone(bset, cone);
2982 static void check_gbr_integer_feasible(struct isl_context_gbr *cgbr)
2984 struct isl_vec *sample;
2986 if (!cgbr->tab)
2987 return;
2989 if (cgbr->tab->empty)
2990 return;
2992 sample = gbr_get_sample(cgbr);
2993 if (!sample)
2994 goto error;
2996 if (sample->size == 0) {
2997 isl_vec_free(sample);
2998 if (isl_tab_mark_empty(cgbr->tab) < 0)
2999 goto error;
3000 return;
3003 if (isl_tab_add_sample(cgbr->tab, sample) < 0)
3004 goto error;
3006 return;
3007 error:
3008 isl_tab_free(cgbr->tab);
3009 cgbr->tab = NULL;
3012 static struct isl_tab *add_gbr_eq(struct isl_tab *tab, isl_int *eq)
3014 if (!tab)
3015 return NULL;
3017 if (isl_tab_extend_cons(tab, 2) < 0)
3018 goto error;
3020 if (isl_tab_add_eq(tab, eq) < 0)
3021 goto error;
3023 return tab;
3024 error:
3025 isl_tab_free(tab);
3026 return NULL;
3029 /* Add the equality described by "eq" to the context.
3030 * If "check" is set, then we check if the context is empty after
3031 * adding the equality.
3032 * If "update" is set, then we check if the samples are still valid.
3034 * We do not explicitly add shifted copies of the equality to
3035 * cgbr->shifted since they would conflict with each other.
3036 * Instead, we directly mark cgbr->shifted empty.
3038 static void context_gbr_add_eq(struct isl_context *context, isl_int *eq,
3039 int check, int update)
3041 struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3043 cgbr->tab = add_gbr_eq(cgbr->tab, eq);
3045 if (cgbr->shifted && !cgbr->shifted->empty && use_shifted(cgbr)) {
3046 if (isl_tab_mark_empty(cgbr->shifted) < 0)
3047 goto error;
3050 if (cgbr->cone && cgbr->cone->n_col != cgbr->cone->n_dead) {
3051 if (isl_tab_extend_cons(cgbr->cone, 2) < 0)
3052 goto error;
3053 if (isl_tab_add_eq(cgbr->cone, eq) < 0)
3054 goto error;
3057 if (check) {
3058 int v = tab_has_valid_sample(cgbr->tab, eq, 1);
3059 if (v < 0)
3060 goto error;
3061 if (!v)
3062 check_gbr_integer_feasible(cgbr);
3064 if (update)
3065 cgbr->tab = check_samples(cgbr->tab, eq, 1);
3066 return;
3067 error:
3068 isl_tab_free(cgbr->tab);
3069 cgbr->tab = NULL;
3072 static void add_gbr_ineq(struct isl_context_gbr *cgbr, isl_int *ineq)
3074 if (!cgbr->tab)
3075 return;
3077 if (isl_tab_extend_cons(cgbr->tab, 1) < 0)
3078 goto error;
3080 if (isl_tab_add_ineq(cgbr->tab, ineq) < 0)
3081 goto error;
3083 if (cgbr->shifted && !cgbr->shifted->empty && use_shifted(cgbr)) {
3084 int i;
3085 unsigned dim;
3086 dim = isl_basic_map_total_dim(cgbr->tab->bmap);
3088 if (isl_tab_extend_cons(cgbr->shifted, 1) < 0)
3089 goto error;
3091 for (i = 0; i < dim; ++i) {
3092 if (!isl_int_is_neg(ineq[1 + i]))
3093 continue;
3094 isl_int_add(ineq[0], ineq[0], ineq[1 + i]);
3097 if (isl_tab_add_ineq(cgbr->shifted, ineq) < 0)
3098 goto error;
3100 for (i = 0; i < dim; ++i) {
3101 if (!isl_int_is_neg(ineq[1 + i]))
3102 continue;
3103 isl_int_sub(ineq[0], ineq[0], ineq[1 + i]);
3107 if (cgbr->cone && cgbr->cone->n_col != cgbr->cone->n_dead) {
3108 if (isl_tab_extend_cons(cgbr->cone, 1) < 0)
3109 goto error;
3110 if (isl_tab_add_ineq(cgbr->cone, ineq) < 0)
3111 goto error;
3114 return;
3115 error:
3116 isl_tab_free(cgbr->tab);
3117 cgbr->tab = NULL;
3120 static void context_gbr_add_ineq(struct isl_context *context, isl_int *ineq,
3121 int check, int update)
3123 struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3125 add_gbr_ineq(cgbr, ineq);
3126 if (!cgbr->tab)
3127 return;
3129 if (check) {
3130 int v = tab_has_valid_sample(cgbr->tab, ineq, 0);
3131 if (v < 0)
3132 goto error;
3133 if (!v)
3134 check_gbr_integer_feasible(cgbr);
3136 if (update)
3137 cgbr->tab = check_samples(cgbr->tab, ineq, 0);
3138 return;
3139 error:
3140 isl_tab_free(cgbr->tab);
3141 cgbr->tab = NULL;
3144 static isl_stat context_gbr_add_ineq_wrap(void *user, isl_int *ineq)
3146 struct isl_context *context = (struct isl_context *)user;
3147 context_gbr_add_ineq(context, ineq, 0, 0);
3148 return context->op->is_ok(context) ? isl_stat_ok : isl_stat_error;
3151 static enum isl_tab_row_sign context_gbr_ineq_sign(struct isl_context *context,
3152 isl_int *ineq, int strict)
3154 struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3155 return tab_ineq_sign(cgbr->tab, ineq, strict);
3158 /* Check whether "ineq" can be added to the tableau without rendering
3159 * it infeasible.
3161 static int context_gbr_test_ineq(struct isl_context *context, isl_int *ineq)
3163 struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3164 struct isl_tab_undo *snap;
3165 struct isl_tab_undo *shifted_snap = NULL;
3166 struct isl_tab_undo *cone_snap = NULL;
3167 int feasible;
3169 if (!cgbr->tab)
3170 return -1;
3172 if (isl_tab_extend_cons(cgbr->tab, 1) < 0)
3173 return -1;
3175 snap = isl_tab_snap(cgbr->tab);
3176 if (cgbr->shifted)
3177 shifted_snap = isl_tab_snap(cgbr->shifted);
3178 if (cgbr->cone)
3179 cone_snap = isl_tab_snap(cgbr->cone);
3180 add_gbr_ineq(cgbr, ineq);
3181 check_gbr_integer_feasible(cgbr);
3182 if (!cgbr->tab)
3183 return -1;
3184 feasible = !cgbr->tab->empty;
3185 if (isl_tab_rollback(cgbr->tab, snap) < 0)
3186 return -1;
3187 if (shifted_snap) {
3188 if (isl_tab_rollback(cgbr->shifted, shifted_snap))
3189 return -1;
3190 } else if (cgbr->shifted) {
3191 isl_tab_free(cgbr->shifted);
3192 cgbr->shifted = NULL;
3194 if (cone_snap) {
3195 if (isl_tab_rollback(cgbr->cone, cone_snap))
3196 return -1;
3197 } else if (cgbr->cone) {
3198 isl_tab_free(cgbr->cone);
3199 cgbr->cone = NULL;
3202 return feasible;
3205 /* Return the column of the last of the variables associated to
3206 * a column that has a non-zero coefficient.
3207 * This function is called in a context where only coefficients
3208 * of parameters or divs can be non-zero.
3210 static int last_non_zero_var_col(struct isl_tab *tab, isl_int *p)
3212 int i;
3213 int col;
3215 if (tab->n_var == 0)
3216 return -1;
3218 for (i = tab->n_var - 1; i >= 0; --i) {
3219 if (i >= tab->n_param && i < tab->n_var - tab->n_div)
3220 continue;
3221 if (tab->var[i].is_row)
3222 continue;
3223 col = tab->var[i].index;
3224 if (!isl_int_is_zero(p[col]))
3225 return col;
3228 return -1;
3231 /* Look through all the recently added equalities in the context
3232 * to see if we can propagate any of them to the main tableau.
3234 * The newly added equalities in the context are encoded as pairs
3235 * of inequalities starting at inequality "first".
3237 * We tentatively add each of these equalities to the main tableau
3238 * and if this happens to result in a row with a final coefficient
3239 * that is one or negative one, we use it to kill a column
3240 * in the main tableau. Otherwise, we discard the tentatively
3241 * added row.
3242 * This tentative addition of equality constraints turns
3243 * on the undo facility of the tableau. Turn it off again
3244 * at the end, assuming it was turned off to begin with.
3246 * Return 0 on success and -1 on failure.
3248 static int propagate_equalities(struct isl_context_gbr *cgbr,
3249 struct isl_tab *tab, unsigned first)
3251 int i;
3252 struct isl_vec *eq = NULL;
3253 isl_bool needs_undo;
3255 needs_undo = isl_tab_need_undo(tab);
3256 if (needs_undo < 0)
3257 goto error;
3258 eq = isl_vec_alloc(tab->mat->ctx, 1 + tab->n_var);
3259 if (!eq)
3260 goto error;
3262 if (isl_tab_extend_cons(tab, (cgbr->tab->bmap->n_ineq - first)/2) < 0)
3263 goto error;
3265 isl_seq_clr(eq->el + 1 + tab->n_param,
3266 tab->n_var - tab->n_param - tab->n_div);
3267 for (i = first; i < cgbr->tab->bmap->n_ineq; i += 2) {
3268 int j;
3269 int r;
3270 struct isl_tab_undo *snap;
3271 snap = isl_tab_snap(tab);
3273 isl_seq_cpy(eq->el, cgbr->tab->bmap->ineq[i], 1 + tab->n_param);
3274 isl_seq_cpy(eq->el + 1 + tab->n_var - tab->n_div,
3275 cgbr->tab->bmap->ineq[i] + 1 + tab->n_param,
3276 tab->n_div);
3278 r = isl_tab_add_row(tab, eq->el);
3279 if (r < 0)
3280 goto error;
3281 r = tab->con[r].index;
3282 j = last_non_zero_var_col(tab, tab->mat->row[r] + 2 + tab->M);
3283 if (j < 0 || j < tab->n_dead ||
3284 !isl_int_is_one(tab->mat->row[r][0]) ||
3285 (!isl_int_is_one(tab->mat->row[r][2 + tab->M + j]) &&
3286 !isl_int_is_negone(tab->mat->row[r][2 + tab->M + j]))) {
3287 if (isl_tab_rollback(tab, snap) < 0)
3288 goto error;
3289 continue;
3291 if (isl_tab_pivot(tab, r, j) < 0)
3292 goto error;
3293 if (isl_tab_kill_col(tab, j) < 0)
3294 goto error;
3296 if (restore_lexmin(tab) < 0)
3297 goto error;
3300 if (!needs_undo)
3301 isl_tab_clear_undo(tab);
3302 isl_vec_free(eq);
3304 return 0;
3305 error:
3306 isl_vec_free(eq);
3307 isl_tab_free(cgbr->tab);
3308 cgbr->tab = NULL;
3309 return -1;
3312 static int context_gbr_detect_equalities(struct isl_context *context,
3313 struct isl_tab *tab)
3315 struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3316 unsigned n_ineq;
3318 if (!cgbr->cone) {
3319 struct isl_basic_set *bset = isl_tab_peek_bset(cgbr->tab);
3320 cgbr->cone = isl_tab_from_recession_cone(bset, 0);
3321 if (!cgbr->cone)
3322 goto error;
3323 if (isl_tab_track_bset(cgbr->cone,
3324 isl_basic_set_copy(bset)) < 0)
3325 goto error;
3327 if (isl_tab_detect_implicit_equalities(cgbr->cone) < 0)
3328 goto error;
3330 n_ineq = cgbr->tab->bmap->n_ineq;
3331 cgbr->tab = isl_tab_detect_equalities(cgbr->tab, cgbr->cone);
3332 if (!cgbr->tab)
3333 return -1;
3334 if (cgbr->tab->bmap->n_ineq > n_ineq &&
3335 propagate_equalities(cgbr, tab, n_ineq) < 0)
3336 return -1;
3338 return 0;
3339 error:
3340 isl_tab_free(cgbr->tab);
3341 cgbr->tab = NULL;
3342 return -1;
3345 static int context_gbr_get_div(struct isl_context *context, struct isl_tab *tab,
3346 struct isl_vec *div)
3348 return get_div(tab, context, div);
3351 static isl_bool context_gbr_insert_div(struct isl_context *context, int pos,
3352 __isl_keep isl_vec *div)
3354 struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3355 if (cgbr->cone) {
3356 int r, n_div, o_div;
3358 n_div = isl_basic_map_dim(cgbr->cone->bmap, isl_dim_div);
3359 o_div = cgbr->cone->n_var - n_div;
3361 if (isl_tab_extend_cons(cgbr->cone, 3) < 0)
3362 return isl_bool_error;
3363 if (isl_tab_extend_vars(cgbr->cone, 1) < 0)
3364 return isl_bool_error;
3365 if ((r = isl_tab_insert_var(cgbr->cone, pos)) <0)
3366 return isl_bool_error;
3368 cgbr->cone->bmap = isl_basic_map_insert_div(cgbr->cone->bmap,
3369 r - o_div, div);
3370 if (!cgbr->cone->bmap)
3371 return isl_bool_error;
3372 if (isl_tab_push_var(cgbr->cone, isl_tab_undo_bmap_div,
3373 &cgbr->cone->var[r]) < 0)
3374 return isl_bool_error;
3376 return context_tab_insert_div(cgbr->tab, pos, div,
3377 context_gbr_add_ineq_wrap, context);
3380 static int context_gbr_best_split(struct isl_context *context,
3381 struct isl_tab *tab)
3383 struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3384 struct isl_tab_undo *snap;
3385 int r;
3387 snap = isl_tab_snap(cgbr->tab);
3388 r = best_split(tab, cgbr->tab);
3390 if (r >= 0 && isl_tab_rollback(cgbr->tab, snap) < 0)
3391 return -1;
3393 return r;
3396 static int context_gbr_is_empty(struct isl_context *context)
3398 struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3399 if (!cgbr->tab)
3400 return -1;
3401 return cgbr->tab->empty;
3404 struct isl_gbr_tab_undo {
3405 struct isl_tab_undo *tab_snap;
3406 struct isl_tab_undo *shifted_snap;
3407 struct isl_tab_undo *cone_snap;
3410 static void *context_gbr_save(struct isl_context *context)
3412 struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3413 struct isl_gbr_tab_undo *snap;
3415 if (!cgbr->tab)
3416 return NULL;
3418 snap = isl_alloc_type(cgbr->tab->mat->ctx, struct isl_gbr_tab_undo);
3419 if (!snap)
3420 return NULL;
3422 snap->tab_snap = isl_tab_snap(cgbr->tab);
3423 if (isl_tab_save_samples(cgbr->tab) < 0)
3424 goto error;
3426 if (cgbr->shifted)
3427 snap->shifted_snap = isl_tab_snap(cgbr->shifted);
3428 else
3429 snap->shifted_snap = NULL;
3431 if (cgbr->cone)
3432 snap->cone_snap = isl_tab_snap(cgbr->cone);
3433 else
3434 snap->cone_snap = NULL;
3436 return snap;
3437 error:
3438 free(snap);
3439 return NULL;
3442 static void context_gbr_restore(struct isl_context *context, void *save)
3444 struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3445 struct isl_gbr_tab_undo *snap = (struct isl_gbr_tab_undo *)save;
3446 if (!snap)
3447 goto error;
3448 if (isl_tab_rollback(cgbr->tab, snap->tab_snap) < 0)
3449 goto error;
3451 if (snap->shifted_snap) {
3452 if (isl_tab_rollback(cgbr->shifted, snap->shifted_snap) < 0)
3453 goto error;
3454 } else if (cgbr->shifted) {
3455 isl_tab_free(cgbr->shifted);
3456 cgbr->shifted = NULL;
3459 if (snap->cone_snap) {
3460 if (isl_tab_rollback(cgbr->cone, snap->cone_snap) < 0)
3461 goto error;
3462 } else if (cgbr->cone) {
3463 isl_tab_free(cgbr->cone);
3464 cgbr->cone = NULL;
3467 free(snap);
3469 return;
3470 error:
3471 free(snap);
3472 isl_tab_free(cgbr->tab);
3473 cgbr->tab = NULL;
3476 static void context_gbr_discard(void *save)
3478 struct isl_gbr_tab_undo *snap = (struct isl_gbr_tab_undo *)save;
3479 free(snap);
3482 static int context_gbr_is_ok(struct isl_context *context)
3484 struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3485 return !!cgbr->tab;
3488 static void context_gbr_invalidate(struct isl_context *context)
3490 struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3491 isl_tab_free(cgbr->tab);
3492 cgbr->tab = NULL;
3495 static __isl_null struct isl_context *context_gbr_free(
3496 struct isl_context *context)
3498 struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3499 isl_tab_free(cgbr->tab);
3500 isl_tab_free(cgbr->shifted);
3501 isl_tab_free(cgbr->cone);
3502 free(cgbr);
3504 return NULL;
3507 struct isl_context_op isl_context_gbr_op = {
3508 context_gbr_detect_nonnegative_parameters,
3509 context_gbr_peek_basic_set,
3510 context_gbr_peek_tab,
3511 context_gbr_add_eq,
3512 context_gbr_add_ineq,
3513 context_gbr_ineq_sign,
3514 context_gbr_test_ineq,
3515 context_gbr_get_div,
3516 context_gbr_insert_div,
3517 context_gbr_detect_equalities,
3518 context_gbr_best_split,
3519 context_gbr_is_empty,
3520 context_gbr_is_ok,
3521 context_gbr_save,
3522 context_gbr_restore,
3523 context_gbr_discard,
3524 context_gbr_invalidate,
3525 context_gbr_free,
3528 static struct isl_context *isl_context_gbr_alloc(__isl_keep isl_basic_set *dom)
3530 struct isl_context_gbr *cgbr;
3532 if (!dom)
3533 return NULL;
3535 cgbr = isl_calloc_type(dom->ctx, struct isl_context_gbr);
3536 if (!cgbr)
3537 return NULL;
3539 cgbr->context.op = &isl_context_gbr_op;
3541 cgbr->shifted = NULL;
3542 cgbr->cone = NULL;
3543 cgbr->tab = isl_tab_from_basic_set(dom, 1);
3544 cgbr->tab = isl_tab_init_samples(cgbr->tab);
3545 if (!cgbr->tab)
3546 goto error;
3547 check_gbr_integer_feasible(cgbr);
3549 return &cgbr->context;
3550 error:
3551 cgbr->context.op->free(&cgbr->context);
3552 return NULL;
3555 /* Allocate a context corresponding to "dom".
3556 * The representation specific fields are initialized by
3557 * isl_context_lex_alloc or isl_context_gbr_alloc.
3558 * The shared "n_unknown" field is initialized to the number
3559 * of final unknown integer divisions in "dom".
3561 static struct isl_context *isl_context_alloc(__isl_keep isl_basic_set *dom)
3563 struct isl_context *context;
3564 int first;
3566 if (!dom)
3567 return NULL;
3569 if (dom->ctx->opt->context == ISL_CONTEXT_LEXMIN)
3570 context = isl_context_lex_alloc(dom);
3571 else
3572 context = isl_context_gbr_alloc(dom);
3574 if (!context)
3575 return NULL;
3577 first = isl_basic_set_first_unknown_div(dom);
3578 if (first < 0)
3579 return context->op->free(context);
3580 context->n_unknown = isl_basic_set_dim(dom, isl_dim_div) - first;
3582 return context;
3585 /* Initialize some common fields of "sol", which keeps track
3586 * of the solution of an optimization problem on "bmap" over
3587 * the domain "dom".
3588 * If "max" is set, then a maximization problem is being solved, rather than
3589 * a minimization problem, which means that the variables in the
3590 * tableau have value "M - x" rather than "M + x".
3592 static isl_stat sol_init(struct isl_sol *sol, __isl_keep isl_basic_map *bmap,
3593 __isl_keep isl_basic_set *dom, int max)
3595 sol->rational = ISL_F_ISSET(bmap, ISL_BASIC_MAP_RATIONAL);
3596 sol->dec_level.callback.run = &sol_dec_level_wrap;
3597 sol->dec_level.sol = sol;
3598 sol->max = max;
3599 sol->n_out = isl_basic_map_dim(bmap, isl_dim_out);
3600 sol->space = isl_basic_map_get_space(bmap);
3602 sol->context = isl_context_alloc(dom);
3603 if (!sol->space || !sol->context)
3604 return isl_stat_error;
3606 return isl_stat_ok;
3609 /* Construct an isl_sol_map structure for accumulating the solution.
3610 * If track_empty is set, then we also keep track of the parts
3611 * of the context where there is no solution.
3612 * If max is set, then we are solving a maximization, rather than
3613 * a minimization problem, which means that the variables in the
3614 * tableau have value "M - x" rather than "M + x".
3616 static struct isl_sol *sol_map_init(__isl_keep isl_basic_map *bmap,
3617 __isl_take isl_basic_set *dom, int track_empty, int max)
3619 struct isl_sol_map *sol_map = NULL;
3620 isl_space *space;
3622 if (!bmap)
3623 goto error;
3625 sol_map = isl_calloc_type(bmap->ctx, struct isl_sol_map);
3626 if (!sol_map)
3627 goto error;
3629 sol_map->sol.free = &sol_map_free;
3630 if (sol_init(&sol_map->sol, bmap, dom, max) < 0)
3631 goto error;
3632 sol_map->sol.add = &sol_map_add_wrap;
3633 sol_map->sol.add_empty = track_empty ? &sol_map_add_empty_wrap : NULL;
3634 space = isl_space_copy(sol_map->sol.space);
3635 sol_map->map = isl_map_alloc_space(space, 1, ISL_MAP_DISJOINT);
3636 if (!sol_map->map)
3637 goto error;
3639 if (track_empty) {
3640 sol_map->empty = isl_set_alloc_space(isl_basic_set_get_space(dom),
3641 1, ISL_SET_DISJOINT);
3642 if (!sol_map->empty)
3643 goto error;
3646 isl_basic_set_free(dom);
3647 return &sol_map->sol;
3648 error:
3649 isl_basic_set_free(dom);
3650 sol_free(&sol_map->sol);
3651 return NULL;
3654 /* Check whether all coefficients of (non-parameter) variables
3655 * are non-positive, meaning that no pivots can be performed on the row.
3657 static int is_critical(struct isl_tab *tab, int row)
3659 int j;
3660 unsigned off = 2 + tab->M;
3662 for (j = tab->n_dead; j < tab->n_col; ++j) {
3663 if (tab->col_var[j] >= 0 &&
3664 (tab->col_var[j] < tab->n_param ||
3665 tab->col_var[j] >= tab->n_var - tab->n_div))
3666 continue;
3668 if (isl_int_is_pos(tab->mat->row[row][off + j]))
3669 return 0;
3672 return 1;
3675 /* Check whether the inequality represented by vec is strict over the integers,
3676 * i.e., there are no integer values satisfying the constraint with
3677 * equality. This happens if the gcd of the coefficients is not a divisor
3678 * of the constant term. If so, scale the constraint down by the gcd
3679 * of the coefficients.
3681 static int is_strict(struct isl_vec *vec)
3683 isl_int gcd;
3684 int strict = 0;
3686 isl_int_init(gcd);
3687 isl_seq_gcd(vec->el + 1, vec->size - 1, &gcd);
3688 if (!isl_int_is_one(gcd)) {
3689 strict = !isl_int_is_divisible_by(vec->el[0], gcd);
3690 isl_int_fdiv_q(vec->el[0], vec->el[0], gcd);
3691 isl_seq_scale_down(vec->el + 1, vec->el + 1, gcd, vec->size-1);
3693 isl_int_clear(gcd);
3695 return strict;
3698 /* Determine the sign of the given row of the main tableau.
3699 * The result is one of
3700 * isl_tab_row_pos: always non-negative; no pivot needed
3701 * isl_tab_row_neg: always non-positive; pivot
3702 * isl_tab_row_any: can be both positive and negative; split
3704 * We first handle some simple cases
3705 * - the row sign may be known already
3706 * - the row may be obviously non-negative
3707 * - the parametric constant may be equal to that of another row
3708 * for which we know the sign. This sign will be either "pos" or
3709 * "any". If it had been "neg" then we would have pivoted before.
3711 * If none of these cases hold, we check the value of the row for each
3712 * of the currently active samples. Based on the signs of these values
3713 * we make an initial determination of the sign of the row.
3715 * all zero -> unk(nown)
3716 * all non-negative -> pos
3717 * all non-positive -> neg
3718 * both negative and positive -> all
3720 * If we end up with "all", we are done.
3721 * Otherwise, we perform a check for positive and/or negative
3722 * values as follows.
3724 * samples neg unk pos
3725 * <0 ? Y N Y N
3726 * pos any pos
3727 * >0 ? Y N Y N
3728 * any neg any neg
3730 * There is no special sign for "zero", because we can usually treat zero
3731 * as either non-negative or non-positive, whatever works out best.
3732 * However, if the row is "critical", meaning that pivoting is impossible
3733 * then we don't want to limp zero with the non-positive case, because
3734 * then we we would lose the solution for those values of the parameters
3735 * where the value of the row is zero. Instead, we treat 0 as non-negative
3736 * ensuring a split if the row can attain both zero and negative values.
3737 * The same happens when the original constraint was one that could not
3738 * be satisfied with equality by any integer values of the parameters.
3739 * In this case, we normalize the constraint, but then a value of zero
3740 * for the normalized constraint is actually a positive value for the
3741 * original constraint, so again we need to treat zero as non-negative.
3742 * In both these cases, we have the following decision tree instead:
3744 * all non-negative -> pos
3745 * all negative -> neg
3746 * both negative and non-negative -> all
3748 * samples neg pos
3749 * <0 ? Y N
3750 * any pos
3751 * >=0 ? Y N
3752 * any neg
3754 static enum isl_tab_row_sign row_sign(struct isl_tab *tab,
3755 struct isl_sol *sol, int row)
3757 struct isl_vec *ineq = NULL;
3758 enum isl_tab_row_sign res = isl_tab_row_unknown;
3759 int critical;
3760 int strict;
3761 int row2;
3763 if (tab->row_sign[row] != isl_tab_row_unknown)
3764 return tab->row_sign[row];
3765 if (is_obviously_nonneg(tab, row))
3766 return isl_tab_row_pos;
3767 for (row2 = tab->n_redundant; row2 < tab->n_row; ++row2) {
3768 if (tab->row_sign[row2] == isl_tab_row_unknown)
3769 continue;
3770 if (identical_parameter_line(tab, row, row2))
3771 return tab->row_sign[row2];
3774 critical = is_critical(tab, row);
3776 ineq = get_row_parameter_ineq(tab, row);
3777 if (!ineq)
3778 goto error;
3780 strict = is_strict(ineq);
3782 res = sol->context->op->ineq_sign(sol->context, ineq->el,
3783 critical || strict);
3785 if (res == isl_tab_row_unknown || res == isl_tab_row_pos) {
3786 /* test for negative values */
3787 int feasible;
3788 isl_seq_neg(ineq->el, ineq->el, ineq->size);
3789 isl_int_sub_ui(ineq->el[0], ineq->el[0], 1);
3791 feasible = sol->context->op->test_ineq(sol->context, ineq->el);
3792 if (feasible < 0)
3793 goto error;
3794 if (!feasible)
3795 res = isl_tab_row_pos;
3796 else
3797 res = (res == isl_tab_row_unknown) ? isl_tab_row_neg
3798 : isl_tab_row_any;
3799 if (res == isl_tab_row_neg) {
3800 isl_seq_neg(ineq->el, ineq->el, ineq->size);
3801 isl_int_sub_ui(ineq->el[0], ineq->el[0], 1);
3805 if (res == isl_tab_row_neg) {
3806 /* test for positive values */
3807 int feasible;
3808 if (!critical && !strict)
3809 isl_int_sub_ui(ineq->el[0], ineq->el[0], 1);
3811 feasible = sol->context->op->test_ineq(sol->context, ineq->el);
3812 if (feasible < 0)
3813 goto error;
3814 if (feasible)
3815 res = isl_tab_row_any;
3818 isl_vec_free(ineq);
3819 return res;
3820 error:
3821 isl_vec_free(ineq);
3822 return isl_tab_row_unknown;
3825 static void find_solutions(struct isl_sol *sol, struct isl_tab *tab);
3827 /* Find solutions for values of the parameters that satisfy the given
3828 * inequality.
3830 * We currently take a snapshot of the context tableau that is reset
3831 * when we return from this function, while we make a copy of the main
3832 * tableau, leaving the original main tableau untouched.
3833 * These are fairly arbitrary choices. Making a copy also of the context
3834 * tableau would obviate the need to undo any changes made to it later,
3835 * while taking a snapshot of the main tableau could reduce memory usage.
3836 * If we were to switch to taking a snapshot of the main tableau,
3837 * we would have to keep in mind that we need to save the row signs
3838 * and that we need to do this before saving the current basis
3839 * such that the basis has been restore before we restore the row signs.
3841 static void find_in_pos(struct isl_sol *sol, struct isl_tab *tab, isl_int *ineq)
3843 void *saved;
3845 if (!sol->context)
3846 goto error;
3847 saved = sol->context->op->save(sol->context);
3849 tab = isl_tab_dup(tab);
3850 if (!tab)
3851 goto error;
3853 sol->context->op->add_ineq(sol->context, ineq, 0, 1);
3855 find_solutions(sol, tab);
3857 if (!sol->error)
3858 sol->context->op->restore(sol->context, saved);
3859 else
3860 sol->context->op->discard(saved);
3861 return;
3862 error:
3863 sol->error = 1;
3866 /* Record the absence of solutions for those values of the parameters
3867 * that do not satisfy the given inequality with equality.
3869 static void no_sol_in_strict(struct isl_sol *sol,
3870 struct isl_tab *tab, struct isl_vec *ineq)
3872 int empty;
3873 void *saved;
3875 if (!sol->context || sol->error)
3876 goto error;
3877 saved = sol->context->op->save(sol->context);
3879 isl_int_sub_ui(ineq->el[0], ineq->el[0], 1);
3881 sol->context->op->add_ineq(sol->context, ineq->el, 1, 0);
3882 if (!sol->context)
3883 goto error;
3885 empty = tab->empty;
3886 tab->empty = 1;
3887 sol_add(sol, tab);
3888 tab->empty = empty;
3890 isl_int_add_ui(ineq->el[0], ineq->el[0], 1);
3892 sol->context->op->restore(sol->context, saved);
3893 return;
3894 error:
3895 sol->error = 1;
3898 /* Reset all row variables that are marked to have a sign that may
3899 * be both positive and negative to have an unknown sign.
3901 static void reset_any_to_unknown(struct isl_tab *tab)
3903 int row;
3905 for (row = tab->n_redundant; row < tab->n_row; ++row) {
3906 if (!isl_tab_var_from_row(tab, row)->is_nonneg)
3907 continue;
3908 if (tab->row_sign[row] == isl_tab_row_any)
3909 tab->row_sign[row] = isl_tab_row_unknown;
3913 /* Compute the lexicographic minimum of the set represented by the main
3914 * tableau "tab" within the context "sol->context_tab".
3915 * On entry the sample value of the main tableau is lexicographically
3916 * less than or equal to this lexicographic minimum.
3917 * Pivots are performed until a feasible point is found, which is then
3918 * necessarily equal to the minimum, or until the tableau is found to
3919 * be infeasible. Some pivots may need to be performed for only some
3920 * feasible values of the context tableau. If so, the context tableau
3921 * is split into a part where the pivot is needed and a part where it is not.
3923 * Whenever we enter the main loop, the main tableau is such that no
3924 * "obvious" pivots need to be performed on it, where "obvious" means
3925 * that the given row can be seen to be negative without looking at
3926 * the context tableau. In particular, for non-parametric problems,
3927 * no pivots need to be performed on the main tableau.
3928 * The caller of find_solutions is responsible for making this property
3929 * hold prior to the first iteration of the loop, while restore_lexmin
3930 * is called before every other iteration.
3932 * Inside the main loop, we first examine the signs of the rows of
3933 * the main tableau within the context of the context tableau.
3934 * If we find a row that is always non-positive for all values of
3935 * the parameters satisfying the context tableau and negative for at
3936 * least one value of the parameters, we perform the appropriate pivot
3937 * and start over. An exception is the case where no pivot can be
3938 * performed on the row. In this case, we require that the sign of
3939 * the row is negative for all values of the parameters (rather than just
3940 * non-positive). This special case is handled inside row_sign, which
3941 * will say that the row can have any sign if it determines that it can
3942 * attain both negative and zero values.
3944 * If we can't find a row that always requires a pivot, but we can find
3945 * one or more rows that require a pivot for some values of the parameters
3946 * (i.e., the row can attain both positive and negative signs), then we split
3947 * the context tableau into two parts, one where we force the sign to be
3948 * non-negative and one where we force is to be negative.
3949 * The non-negative part is handled by a recursive call (through find_in_pos).
3950 * Upon returning from this call, we continue with the negative part and
3951 * perform the required pivot.
3953 * If no such rows can be found, all rows are non-negative and we have
3954 * found a (rational) feasible point. If we only wanted a rational point
3955 * then we are done.
3956 * Otherwise, we check if all values of the sample point of the tableau
3957 * are integral for the variables. If so, we have found the minimal
3958 * integral point and we are done.
3959 * If the sample point is not integral, then we need to make a distinction
3960 * based on whether the constant term is non-integral or the coefficients
3961 * of the parameters. Furthermore, in order to decide how to handle
3962 * the non-integrality, we also need to know whether the coefficients
3963 * of the other columns in the tableau are integral. This leads
3964 * to the following table. The first two rows do not correspond
3965 * to a non-integral sample point and are only mentioned for completeness.
3967 * constant parameters other
3969 * int int int |
3970 * int int rat | -> no problem
3972 * rat int int -> fail
3974 * rat int rat -> cut
3976 * int rat rat |
3977 * rat rat rat | -> parametric cut
3979 * int rat int |
3980 * rat rat int | -> split context
3982 * If the parametric constant is completely integral, then there is nothing
3983 * to be done. If the constant term is non-integral, but all the other
3984 * coefficient are integral, then there is nothing that can be done
3985 * and the tableau has no integral solution.
3986 * If, on the other hand, one or more of the other columns have rational
3987 * coefficients, but the parameter coefficients are all integral, then
3988 * we can perform a regular (non-parametric) cut.
3989 * Finally, if there is any parameter coefficient that is non-integral,
3990 * then we need to involve the context tableau. There are two cases here.
3991 * If at least one other column has a rational coefficient, then we
3992 * can perform a parametric cut in the main tableau by adding a new
3993 * integer division in the context tableau.
3994 * If all other columns have integral coefficients, then we need to
3995 * enforce that the rational combination of parameters (c + \sum a_i y_i)/m
3996 * is always integral. We do this by introducing an integer division
3997 * q = floor((c + \sum a_i y_i)/m) and stipulating that its argument should
3998 * always be integral in the context tableau, i.e., m q = c + \sum a_i y_i.
3999 * Since q is expressed in the tableau as
4000 * c + \sum a_i y_i - m q >= 0
4001 * -c - \sum a_i y_i + m q + m - 1 >= 0
4002 * it is sufficient to add the inequality
4003 * -c - \sum a_i y_i + m q >= 0
4004 * In the part of the context where this inequality does not hold, the
4005 * main tableau is marked as being empty.
4007 static void find_solutions(struct isl_sol *sol, struct isl_tab *tab)
4009 struct isl_context *context;
4010 int r;
4012 if (!tab || sol->error)
4013 goto error;
4015 context = sol->context;
4017 if (tab->empty)
4018 goto done;
4019 if (context->op->is_empty(context))
4020 goto done;
4022 for (r = 0; r >= 0 && tab && !tab->empty; r = restore_lexmin(tab)) {
4023 int flags;
4024 int row;
4025 enum isl_tab_row_sign sgn;
4026 int split = -1;
4027 int n_split = 0;
4029 for (row = tab->n_redundant; row < tab->n_row; ++row) {
4030 if (!isl_tab_var_from_row(tab, row)->is_nonneg)
4031 continue;
4032 sgn = row_sign(tab, sol, row);
4033 if (!sgn)
4034 goto error;
4035 tab->row_sign[row] = sgn;
4036 if (sgn == isl_tab_row_any)
4037 n_split++;
4038 if (sgn == isl_tab_row_any && split == -1)
4039 split = row;
4040 if (sgn == isl_tab_row_neg)
4041 break;
4043 if (row < tab->n_row)
4044 continue;
4045 if (split != -1) {
4046 struct isl_vec *ineq;
4047 if (n_split != 1)
4048 split = context->op->best_split(context, tab);
4049 if (split < 0)
4050 goto error;
4051 ineq = get_row_parameter_ineq(tab, split);
4052 if (!ineq)
4053 goto error;
4054 is_strict(ineq);
4055 reset_any_to_unknown(tab);
4056 tab->row_sign[split] = isl_tab_row_pos;
4057 sol_inc_level(sol);
4058 find_in_pos(sol, tab, ineq->el);
4059 tab->row_sign[split] = isl_tab_row_neg;
4060 isl_seq_neg(ineq->el, ineq->el, ineq->size);
4061 isl_int_sub_ui(ineq->el[0], ineq->el[0], 1);
4062 if (!sol->error)
4063 context->op->add_ineq(context, ineq->el, 0, 1);
4064 isl_vec_free(ineq);
4065 if (sol->error)
4066 goto error;
4067 continue;
4069 if (tab->rational)
4070 break;
4071 row = first_non_integer_row(tab, &flags);
4072 if (row < 0)
4073 break;
4074 if (ISL_FL_ISSET(flags, I_PAR)) {
4075 if (ISL_FL_ISSET(flags, I_VAR)) {
4076 if (isl_tab_mark_empty(tab) < 0)
4077 goto error;
4078 break;
4080 row = add_cut(tab, row);
4081 } else if (ISL_FL_ISSET(flags, I_VAR)) {
4082 struct isl_vec *div;
4083 struct isl_vec *ineq;
4084 int d;
4085 div = get_row_split_div(tab, row);
4086 if (!div)
4087 goto error;
4088 d = context->op->get_div(context, tab, div);
4089 isl_vec_free(div);
4090 if (d < 0)
4091 goto error;
4092 ineq = ineq_for_div(context->op->peek_basic_set(context), d);
4093 if (!ineq)
4094 goto error;
4095 sol_inc_level(sol);
4096 no_sol_in_strict(sol, tab, ineq);
4097 isl_seq_neg(ineq->el, ineq->el, ineq->size);
4098 context->op->add_ineq(context, ineq->el, 1, 1);
4099 isl_vec_free(ineq);
4100 if (sol->error || !context->op->is_ok(context))
4101 goto error;
4102 tab = set_row_cst_to_div(tab, row, d);
4103 if (context->op->is_empty(context))
4104 break;
4105 } else
4106 row = add_parametric_cut(tab, row, context);
4107 if (row < 0)
4108 goto error;
4110 if (r < 0)
4111 goto error;
4112 done:
4113 sol_add(sol, tab);
4114 isl_tab_free(tab);
4115 return;
4116 error:
4117 isl_tab_free(tab);
4118 sol->error = 1;
4121 /* Does "sol" contain a pair of partial solutions that could potentially
4122 * be merged?
4124 * We currently only check that "sol" is not in an error state
4125 * and that there are at least two partial solutions of which the final two
4126 * are defined at the same level.
4128 static int sol_has_mergeable_solutions(struct isl_sol *sol)
4130 if (sol->error)
4131 return 0;
4132 if (!sol->partial)
4133 return 0;
4134 if (!sol->partial->next)
4135 return 0;
4136 return sol->partial->level == sol->partial->next->level;
4139 /* Compute the lexicographic minimum of the set represented by the main
4140 * tableau "tab" within the context "sol->context_tab".
4142 * As a preprocessing step, we first transfer all the purely parametric
4143 * equalities from the main tableau to the context tableau, i.e.,
4144 * parameters that have been pivoted to a row.
4145 * These equalities are ignored by the main algorithm, because the
4146 * corresponding rows may not be marked as being non-negative.
4147 * In parts of the context where the added equality does not hold,
4148 * the main tableau is marked as being empty.
4150 * Before we embark on the actual computation, we save a copy
4151 * of the context. When we return, we check if there are any
4152 * partial solutions that can potentially be merged. If so,
4153 * we perform a rollback to the initial state of the context.
4154 * The merging of partial solutions happens inside calls to
4155 * sol_dec_level that are pushed onto the undo stack of the context.
4156 * If there are no partial solutions that can potentially be merged
4157 * then the rollback is skipped as it would just be wasted effort.
4159 static void find_solutions_main(struct isl_sol *sol, struct isl_tab *tab)
4161 int row;
4162 void *saved;
4164 if (!tab)
4165 goto error;
4167 sol->level = 0;
4169 for (row = tab->n_redundant; row < tab->n_row; ++row) {
4170 int p;
4171 struct isl_vec *eq;
4173 if (tab->row_var[row] < 0)
4174 continue;
4175 if (tab->row_var[row] >= tab->n_param &&
4176 tab->row_var[row] < tab->n_var - tab->n_div)
4177 continue;
4178 if (tab->row_var[row] < tab->n_param)
4179 p = tab->row_var[row];
4180 else
4181 p = tab->row_var[row]
4182 + tab->n_param - (tab->n_var - tab->n_div);
4184 eq = isl_vec_alloc(tab->mat->ctx, 1+tab->n_param+tab->n_div);
4185 if (!eq)
4186 goto error;
4187 get_row_parameter_line(tab, row, eq->el);
4188 isl_int_neg(eq->el[1 + p], tab->mat->row[row][0]);
4189 eq = isl_vec_normalize(eq);
4191 sol_inc_level(sol);
4192 no_sol_in_strict(sol, tab, eq);
4194 isl_seq_neg(eq->el, eq->el, eq->size);
4195 sol_inc_level(sol);
4196 no_sol_in_strict(sol, tab, eq);
4197 isl_seq_neg(eq->el, eq->el, eq->size);
4199 sol->context->op->add_eq(sol->context, eq->el, 1, 1);
4201 isl_vec_free(eq);
4203 if (isl_tab_mark_redundant(tab, row) < 0)
4204 goto error;
4206 if (sol->context->op->is_empty(sol->context))
4207 break;
4209 row = tab->n_redundant - 1;
4212 saved = sol->context->op->save(sol->context);
4214 find_solutions(sol, tab);
4216 if (sol_has_mergeable_solutions(sol))
4217 sol->context->op->restore(sol->context, saved);
4218 else
4219 sol->context->op->discard(saved);
4221 sol->level = 0;
4222 sol_pop(sol);
4224 return;
4225 error:
4226 isl_tab_free(tab);
4227 sol->error = 1;
4230 /* Check if integer division "div" of "dom" also occurs in "bmap".
4231 * If so, return its position within the divs.
4232 * If not, return -1.
4234 static int find_context_div(struct isl_basic_map *bmap,
4235 struct isl_basic_set *dom, unsigned div)
4237 int i;
4238 unsigned b_dim = isl_space_dim(bmap->dim, isl_dim_all);
4239 unsigned d_dim = isl_space_dim(dom->dim, isl_dim_all);
4241 if (isl_int_is_zero(dom->div[div][0]))
4242 return -1;
4243 if (isl_seq_first_non_zero(dom->div[div] + 2 + d_dim, dom->n_div) != -1)
4244 return -1;
4246 for (i = 0; i < bmap->n_div; ++i) {
4247 if (isl_int_is_zero(bmap->div[i][0]))
4248 continue;
4249 if (isl_seq_first_non_zero(bmap->div[i] + 2 + d_dim,
4250 (b_dim - d_dim) + bmap->n_div) != -1)
4251 continue;
4252 if (isl_seq_eq(bmap->div[i], dom->div[div], 2 + d_dim))
4253 return i;
4255 return -1;
4258 /* The correspondence between the variables in the main tableau,
4259 * the context tableau, and the input map and domain is as follows.
4260 * The first n_param and the last n_div variables of the main tableau
4261 * form the variables of the context tableau.
4262 * In the basic map, these n_param variables correspond to the
4263 * parameters and the input dimensions. In the domain, they correspond
4264 * to the parameters and the set dimensions.
4265 * The n_div variables correspond to the integer divisions in the domain.
4266 * To ensure that everything lines up, we may need to copy some of the
4267 * integer divisions of the domain to the map. These have to be placed
4268 * in the same order as those in the context and they have to be placed
4269 * after any other integer divisions that the map may have.
4270 * This function performs the required reordering.
4272 static __isl_give isl_basic_map *align_context_divs(
4273 __isl_take isl_basic_map *bmap, __isl_keep isl_basic_set *dom)
4275 int i;
4276 int common = 0;
4277 int other;
4279 for (i = 0; i < dom->n_div; ++i)
4280 if (find_context_div(bmap, dom, i) != -1)
4281 common++;
4282 other = bmap->n_div - common;
4283 if (dom->n_div - common > 0) {
4284 bmap = isl_basic_map_extend_space(bmap, isl_space_copy(bmap->dim),
4285 dom->n_div - common, 0, 0);
4286 if (!bmap)
4287 return NULL;
4289 for (i = 0; i < dom->n_div; ++i) {
4290 int pos = find_context_div(bmap, dom, i);
4291 if (pos < 0) {
4292 pos = isl_basic_map_alloc_div(bmap);
4293 if (pos < 0)
4294 goto error;
4295 isl_int_set_si(bmap->div[pos][0], 0);
4297 if (pos != other + i)
4298 isl_basic_map_swap_div(bmap, pos, other + i);
4300 return bmap;
4301 error:
4302 isl_basic_map_free(bmap);
4303 return NULL;
4306 /* Base case of isl_tab_basic_map_partial_lexopt, after removing
4307 * some obvious symmetries.
4309 * We make sure the divs in the domain are properly ordered,
4310 * because they will be added one by one in the given order
4311 * during the construction of the solution map.
4312 * Furthermore, make sure that the known integer divisions
4313 * appear before any unknown integer division because the solution
4314 * may depend on the known integer divisions, while anything that
4315 * depends on any variable starting from the first unknown integer
4316 * division is ignored in sol_pma_add.
4318 static struct isl_sol *basic_map_partial_lexopt_base_sol(
4319 __isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom,
4320 __isl_give isl_set **empty, int max,
4321 struct isl_sol *(*init)(__isl_keep isl_basic_map *bmap,
4322 __isl_take isl_basic_set *dom, int track_empty, int max))
4324 struct isl_tab *tab;
4325 struct isl_sol *sol = NULL;
4326 struct isl_context *context;
4328 if (dom->n_div) {
4329 dom = isl_basic_set_sort_divs(dom);
4330 bmap = align_context_divs(bmap, dom);
4332 sol = init(bmap, dom, !!empty, max);
4333 if (!sol)
4334 goto error;
4336 context = sol->context;
4337 if (isl_basic_set_plain_is_empty(context->op->peek_basic_set(context)))
4338 /* nothing */;
4339 else if (isl_basic_map_plain_is_empty(bmap)) {
4340 if (sol->add_empty)
4341 sol->add_empty(sol,
4342 isl_basic_set_copy(context->op->peek_basic_set(context)));
4343 } else {
4344 tab = tab_for_lexmin(bmap,
4345 context->op->peek_basic_set(context), 1, max);
4346 tab = context->op->detect_nonnegative_parameters(context, tab);
4347 find_solutions_main(sol, tab);
4349 if (sol->error)
4350 goto error;
4352 isl_basic_map_free(bmap);
4353 return sol;
4354 error:
4355 sol_free(sol);
4356 isl_basic_map_free(bmap);
4357 return NULL;
4360 /* Base case of isl_tab_basic_map_partial_lexopt, after removing
4361 * some obvious symmetries.
4363 * We call basic_map_partial_lexopt_base_sol and extract the results.
4365 static __isl_give isl_map *basic_map_partial_lexopt_base(
4366 __isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom,
4367 __isl_give isl_set **empty, int max)
4369 isl_map *result = NULL;
4370 struct isl_sol *sol;
4371 struct isl_sol_map *sol_map;
4373 sol = basic_map_partial_lexopt_base_sol(bmap, dom, empty, max,
4374 &sol_map_init);
4375 if (!sol)
4376 return NULL;
4377 sol_map = (struct isl_sol_map *) sol;
4379 result = isl_map_copy(sol_map->map);
4380 if (empty)
4381 *empty = isl_set_copy(sol_map->empty);
4382 sol_free(&sol_map->sol);
4383 return result;
4386 /* Return a count of the number of occurrences of the "n" first
4387 * variables in the inequality constraints of "bmap".
4389 static __isl_give int *count_occurrences(__isl_keep isl_basic_map *bmap,
4390 int n)
4392 int i, j;
4393 isl_ctx *ctx;
4394 int *occurrences;
4396 if (!bmap)
4397 return NULL;
4398 ctx = isl_basic_map_get_ctx(bmap);
4399 occurrences = isl_calloc_array(ctx, int, n);
4400 if (!occurrences)
4401 return NULL;
4403 for (i = 0; i < bmap->n_ineq; ++i) {
4404 for (j = 0; j < n; ++j) {
4405 if (!isl_int_is_zero(bmap->ineq[i][1 + j]))
4406 occurrences[j]++;
4410 return occurrences;
4413 /* Do all of the "n" variables with non-zero coefficients in "c"
4414 * occur in exactly a single constraint.
4415 * "occurrences" is an array of length "n" containing the number
4416 * of occurrences of each of the variables in the inequality constraints.
4418 static int single_occurrence(int n, isl_int *c, int *occurrences)
4420 int i;
4422 for (i = 0; i < n; ++i) {
4423 if (isl_int_is_zero(c[i]))
4424 continue;
4425 if (occurrences[i] != 1)
4426 return 0;
4429 return 1;
4432 /* Do all of the "n" initial variables that occur in inequality constraint
4433 * "ineq" of "bmap" only occur in that constraint?
4435 static int all_single_occurrence(__isl_keep isl_basic_map *bmap, int ineq,
4436 int n)
4438 int i, j;
4440 for (i = 0; i < n; ++i) {
4441 if (isl_int_is_zero(bmap->ineq[ineq][1 + i]))
4442 continue;
4443 for (j = 0; j < bmap->n_ineq; ++j) {
4444 if (j == ineq)
4445 continue;
4446 if (!isl_int_is_zero(bmap->ineq[j][1 + i]))
4447 return 0;
4451 return 1;
4454 /* Structure used during detection of parallel constraints.
4455 * n_in: number of "input" variables: isl_dim_param + isl_dim_in
4456 * n_out: number of "output" variables: isl_dim_out + isl_dim_div
4457 * val: the coefficients of the output variables
4459 struct isl_constraint_equal_info {
4460 unsigned n_in;
4461 unsigned n_out;
4462 isl_int *val;
4465 /* Check whether the coefficients of the output variables
4466 * of the constraint in "entry" are equal to info->val.
4468 static int constraint_equal(const void *entry, const void *val)
4470 isl_int **row = (isl_int **)entry;
4471 const struct isl_constraint_equal_info *info = val;
4473 return isl_seq_eq((*row) + 1 + info->n_in, info->val, info->n_out);
4476 /* Check whether "bmap" has a pair of constraints that have
4477 * the same coefficients for the output variables.
4478 * Note that the coefficients of the existentially quantified
4479 * variables need to be zero since the existentially quantified
4480 * of the result are usually not the same as those of the input.
4481 * Furthermore, check that each of the input variables that occur
4482 * in those constraints does not occur in any other constraint.
4483 * If so, return true and return the row indices of the two constraints
4484 * in *first and *second.
4486 static isl_bool parallel_constraints(__isl_keep isl_basic_map *bmap,
4487 int *first, int *second)
4489 int i;
4490 isl_ctx *ctx;
4491 int *occurrences = NULL;
4492 struct isl_hash_table *table = NULL;
4493 struct isl_hash_table_entry *entry;
4494 struct isl_constraint_equal_info info;
4495 unsigned n_out;
4496 unsigned n_div;
4498 ctx = isl_basic_map_get_ctx(bmap);
4499 table = isl_hash_table_alloc(ctx, bmap->n_ineq);
4500 if (!table)
4501 goto error;
4503 info.n_in = isl_basic_map_dim(bmap, isl_dim_param) +
4504 isl_basic_map_dim(bmap, isl_dim_in);
4505 occurrences = count_occurrences(bmap, info.n_in);
4506 if (info.n_in && !occurrences)
4507 goto error;
4508 n_out = isl_basic_map_dim(bmap, isl_dim_out);
4509 n_div = isl_basic_map_dim(bmap, isl_dim_div);
4510 info.n_out = n_out + n_div;
4511 for (i = 0; i < bmap->n_ineq; ++i) {
4512 uint32_t hash;
4514 info.val = bmap->ineq[i] + 1 + info.n_in;
4515 if (isl_seq_first_non_zero(info.val, n_out) < 0)
4516 continue;
4517 if (isl_seq_first_non_zero(info.val + n_out, n_div) >= 0)
4518 continue;
4519 if (!single_occurrence(info.n_in, bmap->ineq[i] + 1,
4520 occurrences))
4521 continue;
4522 hash = isl_seq_get_hash(info.val, info.n_out);
4523 entry = isl_hash_table_find(ctx, table, hash,
4524 constraint_equal, &info, 1);
4525 if (!entry)
4526 goto error;
4527 if (entry->data)
4528 break;
4529 entry->data = &bmap->ineq[i];
4532 if (i < bmap->n_ineq) {
4533 *first = ((isl_int **)entry->data) - bmap->ineq;
4534 *second = i;
4537 isl_hash_table_free(ctx, table);
4538 free(occurrences);
4540 return i < bmap->n_ineq;
4541 error:
4542 isl_hash_table_free(ctx, table);
4543 free(occurrences);
4544 return isl_bool_error;
4547 /* Given a set of upper bounds in "var", add constraints to "bset"
4548 * that make the i-th bound smallest.
4550 * In particular, if there are n bounds b_i, then add the constraints
4552 * b_i <= b_j for j > i
4553 * b_i < b_j for j < i
4555 static __isl_give isl_basic_set *select_minimum(__isl_take isl_basic_set *bset,
4556 __isl_keep isl_mat *var, int i)
4558 isl_ctx *ctx;
4559 int j, k;
4561 ctx = isl_mat_get_ctx(var);
4563 for (j = 0; j < var->n_row; ++j) {
4564 if (j == i)
4565 continue;
4566 k = isl_basic_set_alloc_inequality(bset);
4567 if (k < 0)
4568 goto error;
4569 isl_seq_combine(bset->ineq[k], ctx->one, var->row[j],
4570 ctx->negone, var->row[i], var->n_col);
4571 isl_int_set_si(bset->ineq[k][var->n_col], 0);
4572 if (j < i)
4573 isl_int_sub_ui(bset->ineq[k][0], bset->ineq[k][0], 1);
4576 bset = isl_basic_set_finalize(bset);
4578 return bset;
4579 error:
4580 isl_basic_set_free(bset);
4581 return NULL;
4584 /* Given a set of upper bounds on the last "input" variable m,
4585 * construct a set that assigns the minimal upper bound to m, i.e.,
4586 * construct a set that divides the space into cells where one
4587 * of the upper bounds is smaller than all the others and assign
4588 * this upper bound to m.
4590 * In particular, if there are n bounds b_i, then the result
4591 * consists of n basic sets, each one of the form
4593 * m = b_i
4594 * b_i <= b_j for j > i
4595 * b_i < b_j for j < i
4597 static __isl_give isl_set *set_minimum(__isl_take isl_space *dim,
4598 __isl_take isl_mat *var)
4600 int i, k;
4601 isl_basic_set *bset = NULL;
4602 isl_set *set = NULL;
4604 if (!dim || !var)
4605 goto error;
4607 set = isl_set_alloc_space(isl_space_copy(dim),
4608 var->n_row, ISL_SET_DISJOINT);
4610 for (i = 0; i < var->n_row; ++i) {
4611 bset = isl_basic_set_alloc_space(isl_space_copy(dim), 0,
4612 1, var->n_row - 1);
4613 k = isl_basic_set_alloc_equality(bset);
4614 if (k < 0)
4615 goto error;
4616 isl_seq_cpy(bset->eq[k], var->row[i], var->n_col);
4617 isl_int_set_si(bset->eq[k][var->n_col], -1);
4618 bset = select_minimum(bset, var, i);
4619 set = isl_set_add_basic_set(set, bset);
4622 isl_space_free(dim);
4623 isl_mat_free(var);
4624 return set;
4625 error:
4626 isl_basic_set_free(bset);
4627 isl_set_free(set);
4628 isl_space_free(dim);
4629 isl_mat_free(var);
4630 return NULL;
4633 /* Given that the last input variable of "bmap" represents the minimum
4634 * of the bounds in "cst", check whether we need to split the domain
4635 * based on which bound attains the minimum.
4637 * A split is needed when the minimum appears in an integer division
4638 * or in an equality. Otherwise, it is only needed if it appears in
4639 * an upper bound that is different from the upper bounds on which it
4640 * is defined.
4642 static isl_bool need_split_basic_map(__isl_keep isl_basic_map *bmap,
4643 __isl_keep isl_mat *cst)
4645 int i, j;
4646 unsigned total;
4647 unsigned pos;
4649 pos = cst->n_col - 1;
4650 total = isl_basic_map_dim(bmap, isl_dim_all);
4652 for (i = 0; i < bmap->n_div; ++i)
4653 if (!isl_int_is_zero(bmap->div[i][2 + pos]))
4654 return isl_bool_true;
4656 for (i = 0; i < bmap->n_eq; ++i)
4657 if (!isl_int_is_zero(bmap->eq[i][1 + pos]))
4658 return isl_bool_true;
4660 for (i = 0; i < bmap->n_ineq; ++i) {
4661 if (isl_int_is_nonneg(bmap->ineq[i][1 + pos]))
4662 continue;
4663 if (!isl_int_is_negone(bmap->ineq[i][1 + pos]))
4664 return isl_bool_true;
4665 if (isl_seq_first_non_zero(bmap->ineq[i] + 1 + pos + 1,
4666 total - pos - 1) >= 0)
4667 return isl_bool_true;
4669 for (j = 0; j < cst->n_row; ++j)
4670 if (isl_seq_eq(bmap->ineq[i], cst->row[j], cst->n_col))
4671 break;
4672 if (j >= cst->n_row)
4673 return isl_bool_true;
4676 return isl_bool_false;
4679 /* Given that the last set variable of "bset" represents the minimum
4680 * of the bounds in "cst", check whether we need to split the domain
4681 * based on which bound attains the minimum.
4683 * We simply call need_split_basic_map here. This is safe because
4684 * the position of the minimum is computed from "cst" and not
4685 * from "bmap".
4687 static isl_bool need_split_basic_set(__isl_keep isl_basic_set *bset,
4688 __isl_keep isl_mat *cst)
4690 return need_split_basic_map(bset_to_bmap(bset), cst);
4693 /* Given that the last set variable of "set" represents the minimum
4694 * of the bounds in "cst", check whether we need to split the domain
4695 * based on which bound attains the minimum.
4697 static isl_bool need_split_set(__isl_keep isl_set *set, __isl_keep isl_mat *cst)
4699 int i;
4701 for (i = 0; i < set->n; ++i) {
4702 isl_bool split;
4704 split = need_split_basic_set(set->p[i], cst);
4705 if (split < 0 || split)
4706 return split;
4709 return isl_bool_false;
4712 /* Given a set of which the last set variable is the minimum
4713 * of the bounds in "cst", split each basic set in the set
4714 * in pieces where one of the bounds is (strictly) smaller than the others.
4715 * This subdivision is given in "min_expr".
4716 * The variable is subsequently projected out.
4718 * We only do the split when it is needed.
4719 * For example if the last input variable m = min(a,b) and the only
4720 * constraints in the given basic set are lower bounds on m,
4721 * i.e., l <= m = min(a,b), then we can simply project out m
4722 * to obtain l <= a and l <= b, without having to split on whether
4723 * m is equal to a or b.
4725 static __isl_give isl_set *split(__isl_take isl_set *empty,
4726 __isl_take isl_set *min_expr, __isl_take isl_mat *cst)
4728 int n_in;
4729 int i;
4730 isl_space *dim;
4731 isl_set *res;
4733 if (!empty || !min_expr || !cst)
4734 goto error;
4736 n_in = isl_set_dim(empty, isl_dim_set);
4737 dim = isl_set_get_space(empty);
4738 dim = isl_space_drop_dims(dim, isl_dim_set, n_in - 1, 1);
4739 res = isl_set_empty(dim);
4741 for (i = 0; i < empty->n; ++i) {
4742 isl_bool split;
4743 isl_set *set;
4745 set = isl_set_from_basic_set(isl_basic_set_copy(empty->p[i]));
4746 split = need_split_basic_set(empty->p[i], cst);
4747 if (split < 0)
4748 set = isl_set_free(set);
4749 else if (split)
4750 set = isl_set_intersect(set, isl_set_copy(min_expr));
4751 set = isl_set_remove_dims(set, isl_dim_set, n_in - 1, 1);
4753 res = isl_set_union_disjoint(res, set);
4756 isl_set_free(empty);
4757 isl_set_free(min_expr);
4758 isl_mat_free(cst);
4759 return res;
4760 error:
4761 isl_set_free(empty);
4762 isl_set_free(min_expr);
4763 isl_mat_free(cst);
4764 return NULL;
4767 /* Given a map of which the last input variable is the minimum
4768 * of the bounds in "cst", split each basic set in the set
4769 * in pieces where one of the bounds is (strictly) smaller than the others.
4770 * This subdivision is given in "min_expr".
4771 * The variable is subsequently projected out.
4773 * The implementation is essentially the same as that of "split".
4775 static __isl_give isl_map *split_domain(__isl_take isl_map *opt,
4776 __isl_take isl_set *min_expr, __isl_take isl_mat *cst)
4778 int n_in;
4779 int i;
4780 isl_space *dim;
4781 isl_map *res;
4783 if (!opt || !min_expr || !cst)
4784 goto error;
4786 n_in = isl_map_dim(opt, isl_dim_in);
4787 dim = isl_map_get_space(opt);
4788 dim = isl_space_drop_dims(dim, isl_dim_in, n_in - 1, 1);
4789 res = isl_map_empty(dim);
4791 for (i = 0; i < opt->n; ++i) {
4792 isl_map *map;
4793 isl_bool split;
4795 map = isl_map_from_basic_map(isl_basic_map_copy(opt->p[i]));
4796 split = need_split_basic_map(opt->p[i], cst);
4797 if (split < 0)
4798 map = isl_map_free(map);
4799 else if (split)
4800 map = isl_map_intersect_domain(map,
4801 isl_set_copy(min_expr));
4802 map = isl_map_remove_dims(map, isl_dim_in, n_in - 1, 1);
4804 res = isl_map_union_disjoint(res, map);
4807 isl_map_free(opt);
4808 isl_set_free(min_expr);
4809 isl_mat_free(cst);
4810 return res;
4811 error:
4812 isl_map_free(opt);
4813 isl_set_free(min_expr);
4814 isl_mat_free(cst);
4815 return NULL;
4818 static __isl_give isl_map *basic_map_partial_lexopt(
4819 __isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom,
4820 __isl_give isl_set **empty, int max);
4822 /* This function is called from basic_map_partial_lexopt_symm.
4823 * The last variable of "bmap" and "dom" corresponds to the minimum
4824 * of the bounds in "cst". "map_space" is the space of the original
4825 * input relation (of basic_map_partial_lexopt_symm) and "set_space"
4826 * is the space of the original domain.
4828 * We recursively call basic_map_partial_lexopt and then plug in
4829 * the definition of the minimum in the result.
4831 static __isl_give isl_map *basic_map_partial_lexopt_symm_core(
4832 __isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom,
4833 __isl_give isl_set **empty, int max, __isl_take isl_mat *cst,
4834 __isl_take isl_space *map_space, __isl_take isl_space *set_space)
4836 isl_map *opt;
4837 isl_set *min_expr;
4839 min_expr = set_minimum(isl_basic_set_get_space(dom), isl_mat_copy(cst));
4841 opt = basic_map_partial_lexopt(bmap, dom, empty, max);
4843 if (empty) {
4844 *empty = split(*empty,
4845 isl_set_copy(min_expr), isl_mat_copy(cst));
4846 *empty = isl_set_reset_space(*empty, set_space);
4849 opt = split_domain(opt, min_expr, cst);
4850 opt = isl_map_reset_space(opt, map_space);
4852 return opt;
4855 /* Extract a domain from "bmap" for the purpose of computing
4856 * a lexicographic optimum.
4858 * This function is only called when the caller wants to compute a full
4859 * lexicographic optimum, i.e., without specifying a domain. In this case,
4860 * the caller is not interested in the part of the domain space where
4861 * there is no solution and the domain can be initialized to those constraints
4862 * of "bmap" that only involve the parameters and the input dimensions.
4863 * This relieves the parametric programming engine from detecting those
4864 * inequalities and transferring them to the context. More importantly,
4865 * it ensures that those inequalities are transferred first and not
4866 * intermixed with inequalities that actually split the domain.
4868 * If the caller does not require the absence of existentially quantified
4869 * variables in the result (i.e., if ISL_OPT_QE is not set in "flags"),
4870 * then the actual domain of "bmap" can be used. This ensures that
4871 * the domain does not need to be split at all just to separate out
4872 * pieces of the domain that do not have a solution from piece that do.
4873 * This domain cannot be used in general because it may involve
4874 * (unknown) existentially quantified variables which will then also
4875 * appear in the solution.
4877 static __isl_give isl_basic_set *extract_domain(__isl_keep isl_basic_map *bmap,
4878 unsigned flags)
4880 int n_div;
4881 int n_out;
4883 n_div = isl_basic_map_dim(bmap, isl_dim_div);
4884 n_out = isl_basic_map_dim(bmap, isl_dim_out);
4885 bmap = isl_basic_map_copy(bmap);
4886 if (ISL_FL_ISSET(flags, ISL_OPT_QE)) {
4887 bmap = isl_basic_map_drop_constraints_involving_dims(bmap,
4888 isl_dim_div, 0, n_div);
4889 bmap = isl_basic_map_drop_constraints_involving_dims(bmap,
4890 isl_dim_out, 0, n_out);
4892 return isl_basic_map_domain(bmap);
4895 #undef TYPE
4896 #define TYPE isl_map
4897 #undef SUFFIX
4898 #define SUFFIX
4899 #include "isl_tab_lexopt_templ.c"
4901 struct isl_sol_for {
4902 struct isl_sol sol;
4903 isl_stat (*fn)(__isl_take isl_basic_set *dom,
4904 __isl_take isl_aff_list *list, void *user);
4905 void *user;
4908 static void sol_for_free(struct isl_sol *sol)
4912 /* Add the solution identified by the tableau and the context tableau.
4913 * In particular, "dom" represents the context and "ma" expresses
4914 * the solution on that context.
4916 * See documentation of sol_add for more details.
4918 * Instead of constructing a basic map, this function calls a user
4919 * defined function with the current context as a basic set and
4920 * a list of affine expressions representing the relation between
4921 * the input and output. The space over which the affine expressions
4922 * are defined is the same as that of the domain. The number of
4923 * affine expressions in the list is equal to the number of output variables.
4925 static void sol_for_add(struct isl_sol_for *sol,
4926 __isl_take isl_basic_set *dom, __isl_take isl_multi_aff *ma)
4928 int i, n;
4929 isl_ctx *ctx;
4930 isl_aff *aff;
4931 isl_aff_list *list;
4933 if (sol->sol.error || !dom || !ma)
4934 goto error;
4936 ctx = isl_basic_set_get_ctx(dom);
4937 n = isl_multi_aff_dim(ma, isl_dim_out);
4938 list = isl_aff_list_alloc(ctx, n);
4939 for (i = 0; i < n; ++i) {
4940 aff = isl_multi_aff_get_aff(ma, i);
4941 list = isl_aff_list_add(list, aff);
4944 dom = isl_basic_set_finalize(dom);
4946 if (sol->fn(isl_basic_set_copy(dom), list, sol->user) < 0)
4947 goto error;
4949 isl_basic_set_free(dom);
4950 isl_multi_aff_free(ma);
4951 return;
4952 error:
4953 isl_basic_set_free(dom);
4954 isl_multi_aff_free(ma);
4955 sol->sol.error = 1;
4958 static void sol_for_add_wrap(struct isl_sol *sol,
4959 __isl_take isl_basic_set *dom, __isl_take isl_multi_aff *ma)
4961 sol_for_add((struct isl_sol_for *)sol, dom, ma);
4964 static struct isl_sol_for *sol_for_init(__isl_keep isl_basic_map *bmap, int max,
4965 isl_stat (*fn)(__isl_take isl_basic_set *dom,
4966 __isl_take isl_aff_list *list, void *user),
4967 void *user)
4969 struct isl_sol_for *sol_for = NULL;
4970 isl_space *dom_dim;
4971 struct isl_basic_set *dom = NULL;
4973 sol_for = isl_calloc_type(bmap->ctx, struct isl_sol_for);
4974 if (!sol_for)
4975 goto error;
4977 dom_dim = isl_space_domain(isl_space_copy(bmap->dim));
4978 dom = isl_basic_set_universe(dom_dim);
4980 sol_for->sol.free = &sol_for_free;
4981 if (sol_init(&sol_for->sol, bmap, dom, max) < 0)
4982 goto error;
4983 sol_for->fn = fn;
4984 sol_for->user = user;
4985 sol_for->sol.add = &sol_for_add_wrap;
4986 sol_for->sol.add_empty = NULL;
4988 isl_basic_set_free(dom);
4989 return sol_for;
4990 error:
4991 isl_basic_set_free(dom);
4992 sol_free(&sol_for->sol);
4993 return NULL;
4996 static void sol_for_find_solutions(struct isl_sol_for *sol_for,
4997 struct isl_tab *tab)
4999 find_solutions_main(&sol_for->sol, tab);
5002 isl_stat isl_basic_map_foreach_lexopt(__isl_keep isl_basic_map *bmap, int max,
5003 isl_stat (*fn)(__isl_take isl_basic_set *dom,
5004 __isl_take isl_aff_list *list, void *user),
5005 void *user)
5007 struct isl_sol_for *sol_for = NULL;
5009 bmap = isl_basic_map_copy(bmap);
5010 bmap = isl_basic_map_detect_equalities(bmap);
5011 if (!bmap)
5012 return isl_stat_error;
5014 sol_for = sol_for_init(bmap, max, fn, user);
5015 if (!sol_for)
5016 goto error;
5018 if (isl_basic_map_plain_is_empty(bmap))
5019 /* nothing */;
5020 else {
5021 struct isl_tab *tab;
5022 struct isl_context *context = sol_for->sol.context;
5023 tab = tab_for_lexmin(bmap,
5024 context->op->peek_basic_set(context), 1, max);
5025 tab = context->op->detect_nonnegative_parameters(context, tab);
5026 sol_for_find_solutions(sol_for, tab);
5027 if (sol_for->sol.error)
5028 goto error;
5031 sol_free(&sol_for->sol);
5032 isl_basic_map_free(bmap);
5033 return isl_stat_ok;
5034 error:
5035 sol_free(&sol_for->sol);
5036 isl_basic_map_free(bmap);
5037 return isl_stat_error;
5040 /* Check if the given sequence of len variables starting at pos
5041 * represents a trivial (i.e., zero) solution.
5042 * The variables are assumed to be non-negative and to come in pairs,
5043 * with each pair representing a variable of unrestricted sign.
5044 * The solution is trivial if each such pair in the sequence consists
5045 * of two identical values, meaning that the variable being represented
5046 * has value zero.
5048 static int region_is_trivial(struct isl_tab *tab, int pos, int len)
5050 int i;
5052 if (len == 0)
5053 return 0;
5055 for (i = 0; i < len; i += 2) {
5056 int neg_row;
5057 int pos_row;
5059 neg_row = tab->var[pos + i].is_row ?
5060 tab->var[pos + i].index : -1;
5061 pos_row = tab->var[pos + i + 1].is_row ?
5062 tab->var[pos + i + 1].index : -1;
5064 if ((neg_row < 0 ||
5065 isl_int_is_zero(tab->mat->row[neg_row][1])) &&
5066 (pos_row < 0 ||
5067 isl_int_is_zero(tab->mat->row[pos_row][1])))
5068 continue;
5070 if (neg_row < 0 || pos_row < 0)
5071 return 0;
5072 if (isl_int_ne(tab->mat->row[neg_row][1],
5073 tab->mat->row[pos_row][1]))
5074 return 0;
5077 return 1;
5080 /* Return the index of the first trivial region or -1 if all regions
5081 * are non-trivial.
5083 static int first_trivial_region(struct isl_tab *tab,
5084 int n_region, struct isl_region *region)
5086 int i;
5088 for (i = 0; i < n_region; ++i) {
5089 if (region_is_trivial(tab, region[i].pos, region[i].len))
5090 return i;
5093 return -1;
5096 /* Check if the solution is optimal, i.e., whether the first
5097 * n_op entries are zero.
5099 static int is_optimal(__isl_keep isl_vec *sol, int n_op)
5101 int i;
5103 for (i = 0; i < n_op; ++i)
5104 if (!isl_int_is_zero(sol->el[1 + i]))
5105 return 0;
5106 return 1;
5109 /* Add constraints to "tab" that ensure that any solution is significantly
5110 * better than that represented by "sol". That is, find the first
5111 * relevant (within first n_op) non-zero coefficient and force it (along
5112 * with all previous coefficients) to be zero.
5113 * If the solution is already optimal (all relevant coefficients are zero),
5114 * then just mark the table as empty.
5115 * "n_zero" is the number of coefficients that have been forced zero
5116 * by previous calls to this function at the same level.
5117 * Return the updated number of forced zero coefficients or -1 on error.
5119 * This function assumes that at least 2 * (n_op - n_zero) more rows and
5120 * at least 2 * (n_op - n_zero) more elements in the constraint array
5121 * are available in the tableau.
5123 static int force_better_solution(struct isl_tab *tab,
5124 __isl_keep isl_vec *sol, int n_op, int n_zero)
5126 int i, n;
5127 isl_ctx *ctx;
5128 isl_vec *v = NULL;
5130 if (!sol)
5131 return -1;
5133 for (i = n_zero; i < n_op; ++i)
5134 if (!isl_int_is_zero(sol->el[1 + i]))
5135 break;
5137 if (i == n_op) {
5138 if (isl_tab_mark_empty(tab) < 0)
5139 return -1;
5140 return n_op;
5143 ctx = isl_vec_get_ctx(sol);
5144 v = isl_vec_alloc(ctx, 1 + tab->n_var);
5145 if (!v)
5146 return -1;
5148 n = i + 1;
5149 for (; i >= n_zero; --i) {
5150 v = isl_vec_clr(v);
5151 isl_int_set_si(v->el[1 + i], -1);
5152 if (add_lexmin_eq(tab, v->el) < 0)
5153 goto error;
5156 isl_vec_free(v);
5157 return n;
5158 error:
5159 isl_vec_free(v);
5160 return -1;
5163 /* Local data at each level of the backtracking procedure of
5164 * isl_tab_basic_set_non_trivial_lexmin.
5166 * "n_zero" is the number of initial coordinates that have already
5167 * been forced to be zero at this level.
5169 struct isl_trivial {
5170 int update;
5171 int n_zero;
5172 int region;
5173 int side;
5174 struct isl_tab_undo *snap;
5177 /* Return the lexicographically smallest non-trivial solution of the
5178 * given ILP problem.
5180 * All variables are assumed to be non-negative.
5182 * n_op is the number of initial coordinates to optimize.
5183 * That is, once a solution has been found, we will only continue looking
5184 * for solutions that result in significantly better values for those
5185 * initial coordinates. That is, we only continue looking for solutions
5186 * that increase the number of initial zeros in this sequence.
5188 * A solution is non-trivial, if it is non-trivial on each of the
5189 * specified regions. Each region represents a sequence of pairs
5190 * of variables. A solution is non-trivial on such a region if
5191 * at least one of these pairs consists of different values, i.e.,
5192 * such that the non-negative variable represented by the pair is non-zero.
5194 * Whenever a conflict is encountered, all constraints involved are
5195 * reported to the caller through a call to "conflict".
5197 * We perform a simple branch-and-bound backtracking search.
5198 * Each level in the search represents an initially trivial region
5199 * that is forced to be non-trivial.
5200 * At each level we consider n cases, where n is the length of the region.
5201 * In terms of the n/2 variables of unrestricted signs being encoded by
5202 * the region, we consider the cases
5203 * x_0 >= 1
5204 * x_0 <= -1
5205 * x_0 = 0 and x_1 >= 1
5206 * x_0 = 0 and x_1 <= -1
5207 * x_0 = 0 and x_1 = 0 and x_2 >= 1
5208 * x_0 = 0 and x_1 = 0 and x_2 <= -1
5209 * ...
5210 * The cases are considered in this order, assuming that each pair
5211 * x_i_a x_i_b represents the value x_i_b - x_i_a.
5212 * That is, x_0 >= 1 is enforced by adding the constraint
5213 * x_0_b - x_0_a >= 1
5215 __isl_give isl_vec *isl_tab_basic_set_non_trivial_lexmin(
5216 __isl_take isl_basic_set *bset, int n_op, int n_region,
5217 struct isl_region *region,
5218 int (*conflict)(int con, void *user), void *user)
5220 int i, j;
5221 int r;
5222 isl_ctx *ctx;
5223 isl_vec *v = NULL;
5224 isl_vec *sol = NULL;
5225 struct isl_tab *tab;
5226 struct isl_trivial *triv = NULL;
5227 int level, init;
5229 if (!bset)
5230 return NULL;
5232 ctx = isl_basic_set_get_ctx(bset);
5233 sol = isl_vec_alloc(ctx, 0);
5235 tab = tab_for_lexmin(bset, NULL, 0, 0);
5236 if (!tab)
5237 goto error;
5238 tab->conflict = conflict;
5239 tab->conflict_user = user;
5241 v = isl_vec_alloc(ctx, 1 + tab->n_var);
5242 triv = isl_calloc_array(ctx, struct isl_trivial, n_region);
5243 if (!v || (n_region && !triv))
5244 goto error;
5246 level = 0;
5247 init = 1;
5249 while (level >= 0) {
5250 int side, base;
5252 if (init) {
5253 tab = cut_to_integer_lexmin(tab, CUT_ONE);
5254 if (!tab)
5255 goto error;
5256 if (tab->empty)
5257 goto backtrack;
5258 r = first_trivial_region(tab, n_region, region);
5259 if (r < 0) {
5260 for (i = 0; i < level; ++i)
5261 triv[i].update = 1;
5262 isl_vec_free(sol);
5263 sol = isl_tab_get_sample_value(tab);
5264 if (!sol)
5265 goto error;
5266 if (is_optimal(sol, n_op))
5267 break;
5268 goto backtrack;
5270 if (level >= n_region)
5271 isl_die(ctx, isl_error_internal,
5272 "nesting level too deep", goto error);
5273 if (isl_tab_extend_cons(tab,
5274 2 * region[r].len + 2 * n_op) < 0)
5275 goto error;
5276 triv[level].region = r;
5277 triv[level].side = 0;
5278 triv[level].update = 0;
5279 triv[level].n_zero = 0;
5282 r = triv[level].region;
5283 side = triv[level].side;
5284 base = 2 * (side/2);
5286 if (side >= region[r].len) {
5287 backtrack:
5288 level--;
5289 init = 0;
5290 if (level >= 0)
5291 if (isl_tab_rollback(tab, triv[level].snap) < 0)
5292 goto error;
5293 continue;
5296 if (triv[level].update) {
5297 triv[level].n_zero = force_better_solution(tab, sol,
5298 n_op, triv[level].n_zero);
5299 if (triv[level].n_zero < 0)
5300 goto error;
5301 triv[level].update = 0;
5304 if (side == base && base >= 2) {
5305 for (j = base - 2; j < base; ++j) {
5306 v = isl_vec_clr(v);
5307 isl_int_set_si(v->el[1 + region[r].pos + j], 1);
5308 if (add_lexmin_eq(tab, v->el) < 0)
5309 goto error;
5313 triv[level].snap = isl_tab_snap(tab);
5314 if (isl_tab_push_basis(tab) < 0)
5315 goto error;
5317 v = isl_vec_clr(v);
5318 isl_int_set_si(v->el[0], -1);
5319 isl_int_set_si(v->el[1 + region[r].pos + side], -1);
5320 isl_int_set_si(v->el[1 + region[r].pos + (side ^ 1)], 1);
5321 tab = add_lexmin_ineq(tab, v->el);
5323 triv[level].side++;
5324 level++;
5325 init = 1;
5328 free(triv);
5329 isl_vec_free(v);
5330 isl_tab_free(tab);
5331 isl_basic_set_free(bset);
5333 return sol;
5334 error:
5335 free(triv);
5336 isl_vec_free(v);
5337 isl_tab_free(tab);
5338 isl_basic_set_free(bset);
5339 isl_vec_free(sol);
5340 return NULL;
5343 /* Wrapper for a tableau that is used for computing
5344 * the lexicographically smallest rational point of a non-negative set.
5345 * This point is represented by the sample value of "tab",
5346 * unless "tab" is empty.
5348 struct isl_tab_lexmin {
5349 isl_ctx *ctx;
5350 struct isl_tab *tab;
5353 /* Free "tl" and return NULL.
5355 __isl_null isl_tab_lexmin *isl_tab_lexmin_free(__isl_take isl_tab_lexmin *tl)
5357 if (!tl)
5358 return NULL;
5359 isl_ctx_deref(tl->ctx);
5360 isl_tab_free(tl->tab);
5361 free(tl);
5363 return NULL;
5366 /* Construct an isl_tab_lexmin for computing
5367 * the lexicographically smallest rational point in "bset",
5368 * assuming that all variables are non-negative.
5370 __isl_give isl_tab_lexmin *isl_tab_lexmin_from_basic_set(
5371 __isl_take isl_basic_set *bset)
5373 isl_ctx *ctx;
5374 isl_tab_lexmin *tl;
5376 if (!bset)
5377 return NULL;
5379 ctx = isl_basic_set_get_ctx(bset);
5380 tl = isl_calloc_type(ctx, struct isl_tab_lexmin);
5381 if (!tl)
5382 goto error;
5383 tl->ctx = ctx;
5384 isl_ctx_ref(ctx);
5385 tl->tab = tab_for_lexmin(bset, NULL, 0, 0);
5386 isl_basic_set_free(bset);
5387 if (!tl->tab)
5388 return isl_tab_lexmin_free(tl);
5389 return tl;
5390 error:
5391 isl_basic_set_free(bset);
5392 isl_tab_lexmin_free(tl);
5393 return NULL;
5396 /* Return the dimension of the set represented by "tl".
5398 int isl_tab_lexmin_dim(__isl_keep isl_tab_lexmin *tl)
5400 return tl ? tl->tab->n_var : -1;
5403 /* Add the equality with coefficients "eq" to "tl", updating the optimal
5404 * solution if needed.
5405 * The equality is added as two opposite inequality constraints.
5407 __isl_give isl_tab_lexmin *isl_tab_lexmin_add_eq(__isl_take isl_tab_lexmin *tl,
5408 isl_int *eq)
5410 unsigned n_var;
5412 if (!tl || !eq)
5413 return isl_tab_lexmin_free(tl);
5415 if (isl_tab_extend_cons(tl->tab, 2) < 0)
5416 return isl_tab_lexmin_free(tl);
5417 n_var = tl->tab->n_var;
5418 isl_seq_neg(eq, eq, 1 + n_var);
5419 tl->tab = add_lexmin_ineq(tl->tab, eq);
5420 isl_seq_neg(eq, eq, 1 + n_var);
5421 tl->tab = add_lexmin_ineq(tl->tab, eq);
5423 if (!tl->tab)
5424 return isl_tab_lexmin_free(tl);
5426 return tl;
5429 /* Return the lexicographically smallest rational point in the basic set
5430 * from which "tl" was constructed.
5431 * If the original input was empty, then return a zero-length vector.
5433 __isl_give isl_vec *isl_tab_lexmin_get_solution(__isl_keep isl_tab_lexmin *tl)
5435 if (!tl)
5436 return NULL;
5437 if (tl->tab->empty)
5438 return isl_vec_alloc(tl->ctx, 0);
5439 else
5440 return isl_tab_get_sample_value(tl->tab);
5443 struct isl_sol_pma {
5444 struct isl_sol sol;
5445 isl_pw_multi_aff *pma;
5446 isl_set *empty;
5449 static void sol_pma_free(struct isl_sol *sol)
5451 struct isl_sol_pma *sol_pma = (struct isl_sol_pma *) sol;
5452 isl_pw_multi_aff_free(sol_pma->pma);
5453 isl_set_free(sol_pma->empty);
5456 /* This function is called for parts of the context where there is
5457 * no solution, with "bset" corresponding to the context tableau.
5458 * Simply add the basic set to the set "empty".
5460 static void sol_pma_add_empty(struct isl_sol_pma *sol,
5461 __isl_take isl_basic_set *bset)
5463 if (!bset || !sol->empty)
5464 goto error;
5466 sol->empty = isl_set_grow(sol->empty, 1);
5467 bset = isl_basic_set_simplify(bset);
5468 bset = isl_basic_set_finalize(bset);
5469 sol->empty = isl_set_add_basic_set(sol->empty, bset);
5470 if (!sol->empty)
5471 sol->sol.error = 1;
5472 return;
5473 error:
5474 isl_basic_set_free(bset);
5475 sol->sol.error = 1;
5478 /* Given a basic set "dom" that represents the context and a tuple of
5479 * affine expressions "maff" defined over this domain, construct
5480 * an isl_pw_multi_aff with a single cell corresponding to "dom" and
5481 * the affine expressions in "maff".
5483 static void sol_pma_add(struct isl_sol_pma *sol,
5484 __isl_take isl_basic_set *dom, __isl_take isl_multi_aff *maff)
5486 isl_pw_multi_aff *pma;
5488 dom = isl_basic_set_simplify(dom);
5489 dom = isl_basic_set_finalize(dom);
5490 pma = isl_pw_multi_aff_alloc(isl_set_from_basic_set(dom), maff);
5491 sol->pma = isl_pw_multi_aff_add_disjoint(sol->pma, pma);
5492 if (!sol->pma)
5493 sol->sol.error = 1;
5496 static void sol_pma_add_empty_wrap(struct isl_sol *sol,
5497 __isl_take isl_basic_set *bset)
5499 sol_pma_add_empty((struct isl_sol_pma *)sol, bset);
5502 static void sol_pma_add_wrap(struct isl_sol *sol,
5503 __isl_take isl_basic_set *dom, __isl_take isl_multi_aff *ma)
5505 sol_pma_add((struct isl_sol_pma *)sol, dom, ma);
5508 /* Construct an isl_sol_pma structure for accumulating the solution.
5509 * If track_empty is set, then we also keep track of the parts
5510 * of the context where there is no solution.
5511 * If max is set, then we are solving a maximization, rather than
5512 * a minimization problem, which means that the variables in the
5513 * tableau have value "M - x" rather than "M + x".
5515 static struct isl_sol *sol_pma_init(__isl_keep isl_basic_map *bmap,
5516 __isl_take isl_basic_set *dom, int track_empty, int max)
5518 struct isl_sol_pma *sol_pma = NULL;
5519 isl_space *space;
5521 if (!bmap)
5522 goto error;
5524 sol_pma = isl_calloc_type(bmap->ctx, struct isl_sol_pma);
5525 if (!sol_pma)
5526 goto error;
5528 sol_pma->sol.free = &sol_pma_free;
5529 if (sol_init(&sol_pma->sol, bmap, dom, max) < 0)
5530 goto error;
5531 sol_pma->sol.add = &sol_pma_add_wrap;
5532 sol_pma->sol.add_empty = track_empty ? &sol_pma_add_empty_wrap : NULL;
5533 space = isl_space_copy(sol_pma->sol.space);
5534 sol_pma->pma = isl_pw_multi_aff_empty(space);
5535 if (!sol_pma->pma)
5536 goto error;
5538 if (track_empty) {
5539 sol_pma->empty = isl_set_alloc_space(isl_basic_set_get_space(dom),
5540 1, ISL_SET_DISJOINT);
5541 if (!sol_pma->empty)
5542 goto error;
5545 isl_basic_set_free(dom);
5546 return &sol_pma->sol;
5547 error:
5548 isl_basic_set_free(dom);
5549 sol_free(&sol_pma->sol);
5550 return NULL;
5553 /* Base case of isl_tab_basic_map_partial_lexopt, after removing
5554 * some obvious symmetries.
5556 * We call basic_map_partial_lexopt_base_sol and extract the results.
5558 static __isl_give isl_pw_multi_aff *basic_map_partial_lexopt_base_pw_multi_aff(
5559 __isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom,
5560 __isl_give isl_set **empty, int max)
5562 isl_pw_multi_aff *result = NULL;
5563 struct isl_sol *sol;
5564 struct isl_sol_pma *sol_pma;
5566 sol = basic_map_partial_lexopt_base_sol(bmap, dom, empty, max,
5567 &sol_pma_init);
5568 if (!sol)
5569 return NULL;
5570 sol_pma = (struct isl_sol_pma *) sol;
5572 result = isl_pw_multi_aff_copy(sol_pma->pma);
5573 if (empty)
5574 *empty = isl_set_copy(sol_pma->empty);
5575 sol_free(&sol_pma->sol);
5576 return result;
5579 /* Given that the last input variable of "maff" represents the minimum
5580 * of some bounds, check whether we need to plug in the expression
5581 * of the minimum.
5583 * In particular, check if the last input variable appears in any
5584 * of the expressions in "maff".
5586 static int need_substitution(__isl_keep isl_multi_aff *maff)
5588 int i;
5589 unsigned pos;
5591 pos = isl_multi_aff_dim(maff, isl_dim_in) - 1;
5593 for (i = 0; i < maff->n; ++i)
5594 if (isl_aff_involves_dims(maff->p[i], isl_dim_in, pos, 1))
5595 return 1;
5597 return 0;
5600 /* Given a set of upper bounds on the last "input" variable m,
5601 * construct a piecewise affine expression that selects
5602 * the minimal upper bound to m, i.e.,
5603 * divide the space into cells where one
5604 * of the upper bounds is smaller than all the others and select
5605 * this upper bound on that cell.
5607 * In particular, if there are n bounds b_i, then the result
5608 * consists of n cell, each one of the form
5610 * b_i <= b_j for j > i
5611 * b_i < b_j for j < i
5613 * The affine expression on this cell is
5615 * b_i
5617 static __isl_give isl_pw_aff *set_minimum_pa(__isl_take isl_space *space,
5618 __isl_take isl_mat *var)
5620 int i;
5621 isl_aff *aff = NULL;
5622 isl_basic_set *bset = NULL;
5623 isl_pw_aff *paff = NULL;
5624 isl_space *pw_space;
5625 isl_local_space *ls = NULL;
5627 if (!space || !var)
5628 goto error;
5630 ls = isl_local_space_from_space(isl_space_copy(space));
5631 pw_space = isl_space_copy(space);
5632 pw_space = isl_space_from_domain(pw_space);
5633 pw_space = isl_space_add_dims(pw_space, isl_dim_out, 1);
5634 paff = isl_pw_aff_alloc_size(pw_space, var->n_row);
5636 for (i = 0; i < var->n_row; ++i) {
5637 isl_pw_aff *paff_i;
5639 aff = isl_aff_alloc(isl_local_space_copy(ls));
5640 bset = isl_basic_set_alloc_space(isl_space_copy(space), 0,
5641 0, var->n_row - 1);
5642 if (!aff || !bset)
5643 goto error;
5644 isl_int_set_si(aff->v->el[0], 1);
5645 isl_seq_cpy(aff->v->el + 1, var->row[i], var->n_col);
5646 isl_int_set_si(aff->v->el[1 + var->n_col], 0);
5647 bset = select_minimum(bset, var, i);
5648 paff_i = isl_pw_aff_alloc(isl_set_from_basic_set(bset), aff);
5649 paff = isl_pw_aff_add_disjoint(paff, paff_i);
5652 isl_local_space_free(ls);
5653 isl_space_free(space);
5654 isl_mat_free(var);
5655 return paff;
5656 error:
5657 isl_aff_free(aff);
5658 isl_basic_set_free(bset);
5659 isl_pw_aff_free(paff);
5660 isl_local_space_free(ls);
5661 isl_space_free(space);
5662 isl_mat_free(var);
5663 return NULL;
5666 /* Given a piecewise multi-affine expression of which the last input variable
5667 * is the minimum of the bounds in "cst", plug in the value of the minimum.
5668 * This minimum expression is given in "min_expr_pa".
5669 * The set "min_expr" contains the same information, but in the form of a set.
5670 * The variable is subsequently projected out.
5672 * The implementation is similar to those of "split" and "split_domain".
5673 * If the variable appears in a given expression, then minimum expression
5674 * is plugged in. Otherwise, if the variable appears in the constraints
5675 * and a split is required, then the domain is split. Otherwise, no split
5676 * is performed.
5678 static __isl_give isl_pw_multi_aff *split_domain_pma(
5679 __isl_take isl_pw_multi_aff *opt, __isl_take isl_pw_aff *min_expr_pa,
5680 __isl_take isl_set *min_expr, __isl_take isl_mat *cst)
5682 int n_in;
5683 int i;
5684 isl_space *space;
5685 isl_pw_multi_aff *res;
5687 if (!opt || !min_expr || !cst)
5688 goto error;
5690 n_in = isl_pw_multi_aff_dim(opt, isl_dim_in);
5691 space = isl_pw_multi_aff_get_space(opt);
5692 space = isl_space_drop_dims(space, isl_dim_in, n_in - 1, 1);
5693 res = isl_pw_multi_aff_empty(space);
5695 for (i = 0; i < opt->n; ++i) {
5696 isl_pw_multi_aff *pma;
5698 pma = isl_pw_multi_aff_alloc(isl_set_copy(opt->p[i].set),
5699 isl_multi_aff_copy(opt->p[i].maff));
5700 if (need_substitution(opt->p[i].maff))
5701 pma = isl_pw_multi_aff_substitute(pma,
5702 isl_dim_in, n_in - 1, min_expr_pa);
5703 else {
5704 isl_bool split;
5705 split = need_split_set(opt->p[i].set, cst);
5706 if (split < 0)
5707 pma = isl_pw_multi_aff_free(pma);
5708 else if (split)
5709 pma = isl_pw_multi_aff_intersect_domain(pma,
5710 isl_set_copy(min_expr));
5712 pma = isl_pw_multi_aff_project_out(pma,
5713 isl_dim_in, n_in - 1, 1);
5715 res = isl_pw_multi_aff_add_disjoint(res, pma);
5718 isl_pw_multi_aff_free(opt);
5719 isl_pw_aff_free(min_expr_pa);
5720 isl_set_free(min_expr);
5721 isl_mat_free(cst);
5722 return res;
5723 error:
5724 isl_pw_multi_aff_free(opt);
5725 isl_pw_aff_free(min_expr_pa);
5726 isl_set_free(min_expr);
5727 isl_mat_free(cst);
5728 return NULL;
5731 static __isl_give isl_pw_multi_aff *basic_map_partial_lexopt_pw_multi_aff(
5732 __isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom,
5733 __isl_give isl_set **empty, int max);
5735 /* This function is called from basic_map_partial_lexopt_symm.
5736 * The last variable of "bmap" and "dom" corresponds to the minimum
5737 * of the bounds in "cst". "map_space" is the space of the original
5738 * input relation (of basic_map_partial_lexopt_symm) and "set_space"
5739 * is the space of the original domain.
5741 * We recursively call basic_map_partial_lexopt and then plug in
5742 * the definition of the minimum in the result.
5744 static __isl_give isl_pw_multi_aff *
5745 basic_map_partial_lexopt_symm_core_pw_multi_aff(
5746 __isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom,
5747 __isl_give isl_set **empty, int max, __isl_take isl_mat *cst,
5748 __isl_take isl_space *map_space, __isl_take isl_space *set_space)
5750 isl_pw_multi_aff *opt;
5751 isl_pw_aff *min_expr_pa;
5752 isl_set *min_expr;
5754 min_expr = set_minimum(isl_basic_set_get_space(dom), isl_mat_copy(cst));
5755 min_expr_pa = set_minimum_pa(isl_basic_set_get_space(dom),
5756 isl_mat_copy(cst));
5758 opt = basic_map_partial_lexopt_pw_multi_aff(bmap, dom, empty, max);
5760 if (empty) {
5761 *empty = split(*empty,
5762 isl_set_copy(min_expr), isl_mat_copy(cst));
5763 *empty = isl_set_reset_space(*empty, set_space);
5766 opt = split_domain_pma(opt, min_expr_pa, min_expr, cst);
5767 opt = isl_pw_multi_aff_reset_space(opt, map_space);
5769 return opt;
5772 #undef TYPE
5773 #define TYPE isl_pw_multi_aff
5774 #undef SUFFIX
5775 #define SUFFIX _pw_multi_aff
5776 #include "isl_tab_lexopt_templ.c"