isl 0.27
[isl.git] / isl_tab_pip.c
blobd78df917583ccd2c2550c268fa279af753402e20
1 /*
2 * Copyright 2008-2009 Katholieke Universiteit Leuven
3 * Copyright 2010 INRIA Saclay
4 * Copyright 2016-2017 Sven Verdoolaege
5 * Copyright 2023 Cerebras Systems
7 * Use of this software is governed by the MIT license
9 * Written by Sven Verdoolaege, K.U.Leuven, Departement
10 * Computerwetenschappen, Celestijnenlaan 200A, B-3001 Leuven, Belgium
11 * and INRIA Saclay - Ile-de-France, Parc Club Orsay Universite,
12 * ZAC des vignes, 4 rue Jacques Monod, 91893 Orsay, France
13 * and Cerebras Systems, 1237 E Arques Ave, Sunnyvale, CA, USA
16 #include <isl_ctx_private.h>
17 #include "isl_map_private.h"
18 #include <isl_seq.h>
19 #include "isl_tab.h"
20 #include "isl_sample.h"
21 #include <isl_mat_private.h>
22 #include <isl_vec_private.h>
23 #include <isl_aff_private.h>
24 #include <isl_constraint_private.h>
25 #include <isl_options_private.h>
26 #include <isl_config.h>
28 #include <bset_to_bmap.c>
31 * The implementation of parametric integer linear programming in this file
32 * was inspired by the paper "Parametric Integer Programming" and the
33 * report "Solving systems of affine (in)equalities" by Paul Feautrier
34 * (and others).
36 * The strategy used for obtaining a feasible solution is different
37 * from the one used in isl_tab.c. In particular, in isl_tab.c,
38 * upon finding a constraint that is not yet satisfied, we pivot
39 * in a row that increases the constant term of the row holding the
40 * constraint, making sure the sample solution remains feasible
41 * for all the constraints it already satisfied.
42 * Here, we always pivot in the row holding the constraint,
43 * choosing a column that induces the lexicographically smallest
44 * increment to the sample solution.
46 * By starting out from a sample value that is lexicographically
47 * smaller than any integer point in the problem space, the first
48 * feasible integer sample point we find will also be the lexicographically
49 * smallest. If all variables can be assumed to be non-negative,
50 * then the initial sample value may be chosen equal to zero.
51 * However, we will not make this assumption. Instead, we apply
52 * the "big parameter" trick. Any variable x is then not directly
53 * used in the tableau, but instead it is represented by another
54 * variable x' = M + x, where M is an arbitrarily large (positive)
55 * value. x' is therefore always non-negative, whatever the value of x.
56 * Taking as initial sample value x' = 0 corresponds to x = -M,
57 * which is always smaller than any possible value of x.
59 * The big parameter trick is used in the main tableau and
60 * also in the context tableau if isl_context_lex is used.
61 * In this case, each tableaus has its own big parameter.
62 * Before doing any real work, we check if all the parameters
63 * happen to be non-negative. If so, we drop the column corresponding
64 * to M from the initial context tableau.
65 * If isl_context_gbr is used, then the big parameter trick is only
66 * used in the main tableau.
69 struct isl_context;
70 struct isl_context_op {
71 /* detect nonnegative parameters in context and mark them in tab */
72 struct isl_tab *(*detect_nonnegative_parameters)(
73 struct isl_context *context, struct isl_tab *tab);
74 /* return temporary reference to basic set representation of context */
75 struct isl_basic_set *(*peek_basic_set)(struct isl_context *context);
76 /* return temporary reference to tableau representation of context */
77 struct isl_tab *(*peek_tab)(struct isl_context *context);
78 /* add equality; check is 1 if eq may not be valid;
79 * update is 1 if we may want to call ineq_sign on context later.
81 void (*add_eq)(struct isl_context *context, isl_int *eq,
82 int check, int update);
83 /* add inequality; check is 1 if ineq may not be valid;
84 * update is 1 if we may want to call ineq_sign on context later.
86 void (*add_ineq)(struct isl_context *context, isl_int *ineq,
87 int check, int update);
88 /* check sign of ineq based on previous information.
89 * strict is 1 if saturation should be treated as a positive sign.
91 enum isl_tab_row_sign (*ineq_sign)(struct isl_context *context,
92 isl_int *ineq, int strict);
93 /* check if inequality maintains feasibility */
94 int (*test_ineq)(struct isl_context *context, isl_int *ineq);
95 /* return index of a div that corresponds to "div" */
96 int (*get_div)(struct isl_context *context, struct isl_tab *tab,
97 struct isl_vec *div);
98 /* insert div "div" to context at "pos" and return non-negativity */
99 isl_bool (*insert_div)(struct isl_context *context, int pos,
100 __isl_keep isl_vec *div);
101 int (*detect_equalities)(struct isl_context *context,
102 struct isl_tab *tab);
103 /* return row index of "best" split */
104 int (*best_split)(struct isl_context *context, struct isl_tab *tab);
105 /* check if context has already been determined to be empty */
106 int (*is_empty)(struct isl_context *context);
107 /* check if context is still usable */
108 int (*is_ok)(struct isl_context *context);
109 /* save a copy/snapshot of context */
110 void *(*save)(struct isl_context *context);
111 /* restore saved context */
112 void (*restore)(struct isl_context *context, void *);
113 /* discard saved context */
114 void (*discard)(void *);
115 /* invalidate context */
116 void (*invalidate)(struct isl_context *context);
117 /* free context */
118 __isl_null struct isl_context *(*free)(struct isl_context *context);
121 /* Shared parts of context representation.
123 * "n_unknown" is the number of final unknown integer divisions
124 * in the input domain.
126 struct isl_context {
127 struct isl_context_op *op;
128 int n_unknown;
131 struct isl_context_lex {
132 struct isl_context context;
133 struct isl_tab *tab;
136 /* A stack (linked list) of solutions of subtrees of the search space.
138 * "ma" describes the solution as a function of "dom".
139 * In particular, the domain space of "ma" is equal to the space of "dom".
141 * If "ma" is NULL, then there is no solution on "dom".
143 struct isl_partial_sol {
144 int level;
145 struct isl_basic_set *dom;
146 isl_multi_aff *ma;
148 struct isl_partial_sol *next;
151 struct isl_sol;
152 struct isl_sol_callback {
153 struct isl_tab_callback callback;
154 struct isl_sol *sol;
157 /* isl_sol is an interface for constructing a solution to
158 * a parametric integer linear programming problem.
159 * Every time the algorithm reaches a state where a solution
160 * can be read off from the tableau, the function "add" is called
161 * on the isl_sol passed to find_solutions_main. In a state where
162 * the tableau is empty, "add_empty" is called instead.
163 * "free" is called to free the implementation specific fields, if any.
165 * "error" is set if some error has occurred. This flag invalidates
166 * the remainder of the data structure.
167 * If "rational" is set, then a rational optimization is being performed.
168 * "level" is the current level in the tree with nodes for each
169 * split in the context.
170 * If "max" is set, then a maximization problem is being solved, rather than
171 * a minimization problem, which means that the variables in the
172 * tableau have value "M - x" rather than "M + x".
173 * "n_out" is the number of output dimensions in the input.
174 * "space" is the space in which the solution (and also the input) lives.
176 * The context tableau is owned by isl_sol and is updated incrementally.
178 * There are currently two implementations of this interface,
179 * isl_sol_map, which simply collects the solutions in an isl_map
180 * and (optionally) the parts of the context where there is no solution
181 * in an isl_set, and
182 * isl_sol_pma, which collects an isl_pw_multi_aff instead.
184 struct isl_sol {
185 int error;
186 int rational;
187 int level;
188 int max;
189 isl_size 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 /* Add equality constraint "eq" to the context of "sol".
219 * "check" is set if "eq" is not known to be a valid constraint.
220 * "update" is set if ineq_sign() may still get called on the context.
222 static void sol_context_add_eq(struct isl_sol *sol, isl_int *eq, int check,
223 int update)
225 sol->context->op->add_eq(sol->context, eq, check, update);
226 if (!sol->context->op->is_ok(sol->context))
227 sol->error = 1;
230 /* Add inequality constraint "ineq" to the context of "sol".
231 * "check" is set if "ineq" is not known to be a valid constraint.
232 * "update" is set if ineq_sign() may still get called on the context.
234 static void sol_context_add_ineq(struct isl_sol *sol, isl_int *ineq, int check,
235 int update)
237 if (sol->error)
238 return;
239 sol->context->op->add_ineq(sol->context, ineq, check, update);
240 if (!sol->context->op->is_ok(sol->context))
241 sol->error = 1;
244 /* Push a partial solution represented by a domain and function "ma"
245 * onto the stack of partial solutions.
246 * If "ma" is NULL, then "dom" represents a part of the domain
247 * with no solution.
249 static void sol_push_sol(struct isl_sol *sol,
250 __isl_take isl_basic_set *dom, __isl_take isl_multi_aff *ma)
252 struct isl_partial_sol *partial;
254 if (sol->error || !dom)
255 goto error;
257 partial = isl_alloc_type(dom->ctx, struct isl_partial_sol);
258 if (!partial)
259 goto error;
261 partial->level = sol->level;
262 partial->dom = dom;
263 partial->ma = ma;
264 partial->next = sol->partial;
266 sol->partial = partial;
268 return;
269 error:
270 isl_basic_set_free(dom);
271 isl_multi_aff_free(ma);
272 sol->error = 1;
275 /* Check that the final columns of "M", starting at "first", are zero.
277 static isl_stat check_final_columns_are_zero(__isl_keep isl_mat *M,
278 unsigned first)
280 int i;
281 isl_size rows, cols;
282 unsigned n;
284 rows = isl_mat_rows(M);
285 cols = isl_mat_cols(M);
286 if (rows < 0 || cols < 0)
287 return isl_stat_error;
288 n = cols - first;
289 for (i = 0; i < rows; ++i)
290 if (isl_seq_any_non_zero(M->row[i] + first, n))
291 isl_die(isl_mat_get_ctx(M), isl_error_internal,
292 "final columns should be zero",
293 return isl_stat_error);
294 return isl_stat_ok;
297 /* Set the affine expressions in "ma" according to the rows in "M", which
298 * are defined over the local space "ls".
299 * The matrix "M" may have extra (zero) columns beyond the number
300 * of variables in "ls".
302 static __isl_give isl_multi_aff *set_from_affine_matrix(
303 __isl_take isl_multi_aff *ma, __isl_take isl_local_space *ls,
304 __isl_take isl_mat *M)
306 int i;
307 isl_size dim;
308 isl_aff *aff;
310 dim = isl_local_space_dim(ls, isl_dim_all);
311 if (!ma || dim < 0 || !M)
312 goto error;
314 if (check_final_columns_are_zero(M, 1 + dim) < 0)
315 goto error;
316 for (i = 1; i < M->n_row; ++i) {
317 aff = isl_aff_alloc(isl_local_space_copy(ls));
318 if (aff) {
319 isl_int_set(aff->v->el[0], M->row[0][0]);
320 isl_seq_cpy(aff->v->el + 1, M->row[i], 1 + dim);
322 aff = isl_aff_normalize(aff);
323 ma = isl_multi_aff_set_aff(ma, i - 1, aff);
325 isl_local_space_free(ls);
326 isl_mat_free(M);
328 return ma;
329 error:
330 isl_local_space_free(ls);
331 isl_mat_free(M);
332 isl_multi_aff_free(ma);
333 return NULL;
336 /* Push a partial solution represented by a domain and mapping M
337 * onto the stack of partial solutions.
339 * The affine matrix "M" maps the dimensions of the context
340 * to the output variables. Convert it into an isl_multi_aff and
341 * then call sol_push_sol.
343 * Note that the description of the initial context may have involved
344 * existentially quantified variables, in which case they also appear
345 * in "dom". These need to be removed before creating the affine
346 * expression because an affine expression cannot be defined in terms
347 * of existentially quantified variables without a known representation.
348 * Since newly added integer divisions are inserted before these
349 * existentially quantified variables, they are still in the final
350 * positions and the corresponding final columns of "M" are zero
351 * because align_context_divs adds the existentially quantified
352 * variables of the context to the main tableau without any constraints and
353 * any equality constraints that are added later on can only serve
354 * to eliminate these existentially quantified variables.
356 static void sol_push_sol_mat(struct isl_sol *sol,
357 __isl_take isl_basic_set *dom, __isl_take isl_mat *M)
359 isl_local_space *ls;
360 isl_multi_aff *ma;
361 isl_size n_div;
362 int n_known;
364 n_div = isl_basic_set_dim(dom, isl_dim_div);
365 if (n_div < 0)
366 goto error;
367 n_known = n_div - sol->context->n_unknown;
369 ma = isl_multi_aff_alloc(isl_space_copy(sol->space));
370 ls = isl_basic_set_get_local_space(dom);
371 ls = isl_local_space_drop_dims(ls, isl_dim_div,
372 n_known, n_div - n_known);
373 ma = set_from_affine_matrix(ma, ls, M);
375 if (!ma)
376 dom = isl_basic_set_free(dom);
377 sol_push_sol(sol, dom, ma);
378 return;
379 error:
380 isl_basic_set_free(dom);
381 isl_mat_free(M);
382 sol_push_sol(sol, NULL, NULL);
385 /* Pop one partial solution from the partial solution stack and
386 * pass it on to sol->add or sol->add_empty.
388 static void sol_pop_one(struct isl_sol *sol)
390 struct isl_partial_sol *partial;
392 partial = sol->partial;
393 sol->partial = partial->next;
395 if (partial->ma)
396 sol->add(sol, partial->dom, partial->ma);
397 else
398 sol->add_empty(sol, partial->dom);
399 free(partial);
402 /* Return a fresh copy of the domain represented by the context tableau.
404 static struct isl_basic_set *sol_domain(struct isl_sol *sol)
406 struct isl_basic_set *bset;
408 if (sol->error)
409 return NULL;
411 bset = isl_basic_set_dup(sol->context->op->peek_basic_set(sol->context));
412 bset = isl_basic_set_update_from_tab(bset,
413 sol->context->op->peek_tab(sol->context));
415 return bset;
418 /* Check whether two partial solutions have the same affine expressions.
420 static isl_bool same_solution(struct isl_partial_sol *s1,
421 struct isl_partial_sol *s2)
423 if (!s1->ma != !s2->ma)
424 return isl_bool_false;
425 if (!s1->ma)
426 return isl_bool_true;
428 return isl_multi_aff_plain_is_equal(s1->ma, s2->ma);
431 /* Swap the initial two partial solutions in "sol".
433 * That is, go from
435 * sol->partial = p1; p1->next = p2; p2->next = p3
437 * to
439 * sol->partial = p2; p2->next = p1; p1->next = p3
441 static void swap_initial(struct isl_sol *sol)
443 struct isl_partial_sol *partial;
445 partial = sol->partial;
446 sol->partial = partial->next;
447 partial->next = partial->next->next;
448 sol->partial->next = partial;
451 /* Combine the initial two partial solution of "sol" into
452 * a partial solution with the current context domain of "sol" and
453 * the function description of the second partial solution in the list.
454 * The level of the new partial solution is set to the current level.
456 * That is, the first two partial solutions (D1,M1) and (D2,M2) are
457 * replaced by (D,M2), where D is the domain of "sol", which is assumed
458 * to be the union of D1 and D2, while M1 is assumed to be equal to M2
459 * (at least on D1).
461 static isl_stat combine_initial_into_second(struct isl_sol *sol)
463 struct isl_partial_sol *partial;
464 isl_basic_set *bset;
466 partial = sol->partial;
468 bset = sol_domain(sol);
469 isl_basic_set_free(partial->next->dom);
470 partial->next->dom = bset;
471 partial->next->level = sol->level;
473 if (!bset)
474 return isl_stat_error;
476 sol->partial = partial->next;
477 isl_basic_set_free(partial->dom);
478 isl_multi_aff_free(partial->ma);
479 free(partial);
481 return isl_stat_ok;
484 /* Are "ma1" and "ma2" equal to each other on "dom"?
486 * Combine "ma1" and "ma2" with "dom" and check if the results are the same.
487 * "dom" may have existentially quantified variables. Eliminate them first
488 * as otherwise they would have to be eliminated twice, in a more complicated
489 * context.
491 static isl_bool equal_on_domain(__isl_keep isl_multi_aff *ma1,
492 __isl_keep isl_multi_aff *ma2, __isl_keep isl_basic_set *dom)
494 isl_set *set;
495 isl_pw_multi_aff *pma1, *pma2;
496 isl_bool equal;
498 set = isl_basic_set_compute_divs(isl_basic_set_copy(dom));
499 pma1 = isl_pw_multi_aff_alloc(isl_set_copy(set),
500 isl_multi_aff_copy(ma1));
501 pma2 = isl_pw_multi_aff_alloc(set, isl_multi_aff_copy(ma2));
502 equal = isl_pw_multi_aff_is_equal(pma1, pma2);
503 isl_pw_multi_aff_free(pma1);
504 isl_pw_multi_aff_free(pma2);
506 return equal;
509 /* The initial two partial solutions of "sol" are known to be at
510 * the same level.
511 * If they represent the same solution (on different parts of the domain),
512 * then combine them into a single solution at the current level.
513 * Otherwise, pop them both.
515 * Even if the two partial solution are not obviously the same,
516 * one may still be a simplification of the other over its own domain.
517 * Also check if the two sets of affine functions are equal when
518 * restricted to one of the domains. If so, combine the two
519 * using the set of affine functions on the other domain.
520 * That is, for two partial solutions (D1,M1) and (D2,M2),
521 * if M1 = M2 on D1, then the pair of partial solutions can
522 * be replaced by (D1+D2,M2) and similarly when M1 = M2 on D2.
524 static isl_stat combine_initial_if_equal(struct isl_sol *sol)
526 struct isl_partial_sol *partial;
527 isl_bool same;
529 partial = sol->partial;
531 same = same_solution(partial, partial->next);
532 if (same < 0)
533 return isl_stat_error;
534 if (same)
535 return combine_initial_into_second(sol);
536 if (partial->ma && partial->next->ma) {
537 same = equal_on_domain(partial->ma, partial->next->ma,
538 partial->dom);
539 if (same < 0)
540 return isl_stat_error;
541 if (same)
542 return combine_initial_into_second(sol);
543 same = equal_on_domain(partial->ma, partial->next->ma,
544 partial->next->dom);
545 if (same) {
546 swap_initial(sol);
547 return combine_initial_into_second(sol);
551 sol_pop_one(sol);
552 sol_pop_one(sol);
554 return isl_stat_ok;
557 /* Pop all solutions from the partial solution stack that were pushed onto
558 * the stack at levels that are deeper than the current level.
559 * If the two topmost elements on the stack have the same level
560 * and represent the same solution, then their domains are combined.
561 * This combined domain is the same as the current context domain
562 * as sol_pop is called each time we move back to a higher level.
563 * If the outer level (0) has been reached, then all partial solutions
564 * at the current level are also popped off.
566 static void sol_pop(struct isl_sol *sol)
568 struct isl_partial_sol *partial;
570 if (sol->error)
571 return;
573 partial = sol->partial;
574 if (!partial)
575 return;
577 if (partial->level == 0 && sol->level == 0) {
578 for (partial = sol->partial; partial; partial = sol->partial)
579 sol_pop_one(sol);
580 return;
583 if (partial->level <= sol->level)
584 return;
586 if (partial->next && partial->next->level == partial->level) {
587 if (combine_initial_if_equal(sol) < 0)
588 goto error;
589 } else
590 sol_pop_one(sol);
592 if (sol->level == 0) {
593 for (partial = sol->partial; partial; partial = sol->partial)
594 sol_pop_one(sol);
595 return;
598 if (0)
599 error: sol->error = 1;
602 static void sol_dec_level(struct isl_sol *sol)
604 if (sol->error)
605 return;
607 sol->level--;
609 sol_pop(sol);
612 static isl_stat sol_dec_level_wrap(struct isl_tab_callback *cb)
614 struct isl_sol_callback *callback = (struct isl_sol_callback *)cb;
616 sol_dec_level(callback->sol);
618 return callback->sol->error ? isl_stat_error : isl_stat_ok;
621 /* Move down to next level and push callback onto context tableau
622 * to decrease the level again when it gets rolled back across
623 * the current state. That is, dec_level will be called with
624 * the context tableau in the same state as it is when inc_level
625 * is called.
627 static void sol_inc_level(struct isl_sol *sol)
629 struct isl_tab *tab;
631 if (sol->error)
632 return;
634 sol->level++;
635 tab = sol->context->op->peek_tab(sol->context);
636 if (isl_tab_push_callback(tab, &sol->dec_level.callback) < 0)
637 sol->error = 1;
640 static void scale_rows(struct isl_mat *mat, isl_int m, int n_row)
642 int i;
644 if (isl_int_is_one(m))
645 return;
647 for (i = 0; i < n_row; ++i)
648 isl_seq_scale(mat->row[i], mat->row[i], m, mat->n_col);
651 /* Add the solution identified by the tableau and the context tableau.
653 * The layout of the variables is as follows.
654 * tab->n_var is equal to the total number of variables in the input
655 * map (including divs that were copied from the context)
656 * + the number of extra divs constructed
657 * Of these, the first tab->n_param and the last tab->n_div variables
658 * correspond to the variables in the context, i.e.,
659 * tab->n_param + tab->n_div = context_tab->n_var
660 * tab->n_param is equal to the number of parameters and input
661 * dimensions in the input map
662 * tab->n_div is equal to the number of divs in the context
664 * If there is no solution, then call add_empty with a basic set
665 * that corresponds to the context tableau. (If add_empty is NULL,
666 * then do nothing).
668 * If there is a solution, then first construct a matrix that maps
669 * all dimensions of the context to the output variables, i.e.,
670 * the output dimensions in the input map.
671 * The divs in the input map (if any) that do not correspond to any
672 * div in the context do not appear in the solution.
673 * The algorithm will make sure that they have an integer value,
674 * but these values themselves are of no interest.
675 * We have to be careful not to drop or rearrange any divs in the
676 * context because that would change the meaning of the matrix.
678 * To extract the value of the output variables, it should be noted
679 * that we always use a big parameter M in the main tableau and so
680 * the variable stored in this tableau is not an output variable x itself, but
681 * x' = M + x (in case of minimization)
682 * or
683 * x' = M - x (in case of maximization)
684 * If x' appears in a column, then its optimal value is zero,
685 * which means that the optimal value of x is an unbounded number
686 * (-M for minimization and M for maximization).
687 * We currently assume that the output dimensions in the original map
688 * are bounded, so this cannot occur.
689 * Similarly, when x' appears in a row, then the coefficient of M in that
690 * row is necessarily 1.
691 * If the row in the tableau represents
692 * d x' = c + d M + e(y)
693 * then, in case of minimization, the corresponding row in the matrix
694 * will be
695 * a c + a e(y)
696 * with a d = m, the (updated) common denominator of the matrix.
697 * In case of maximization, the row will be
698 * -a c - a e(y)
700 static void sol_add(struct isl_sol *sol, struct isl_tab *tab)
702 struct isl_basic_set *bset = NULL;
703 struct isl_mat *mat = NULL;
704 unsigned off;
705 int row;
706 isl_int m;
708 if (sol->error || !tab)
709 goto error;
711 if (tab->empty && !sol->add_empty)
712 return;
713 if (sol->context->op->is_empty(sol->context))
714 return;
716 bset = sol_domain(sol);
718 if (tab->empty) {
719 sol_push_sol(sol, bset, NULL);
720 return;
723 off = 2 + tab->M;
725 mat = isl_mat_alloc(tab->mat->ctx, 1 + sol->n_out,
726 1 + tab->n_param + tab->n_div);
727 if (!mat)
728 goto error;
730 isl_int_init(m);
732 isl_seq_clr(mat->row[0] + 1, mat->n_col - 1);
733 isl_int_set_si(mat->row[0][0], 1);
734 for (row = 0; row < sol->n_out; ++row) {
735 int i = tab->n_param + row;
736 int r, j;
738 isl_seq_clr(mat->row[1 + row], mat->n_col);
739 if (!tab->var[i].is_row) {
740 if (tab->M)
741 isl_die(mat->ctx, isl_error_invalid,
742 "unbounded optimum", goto error2);
743 continue;
746 r = tab->var[i].index;
747 if (tab->M &&
748 isl_int_ne(tab->mat->row[r][2], tab->mat->row[r][0]))
749 isl_die(mat->ctx, isl_error_invalid,
750 "unbounded optimum", goto error2);
751 isl_int_gcd(m, mat->row[0][0], tab->mat->row[r][0]);
752 isl_int_divexact(m, tab->mat->row[r][0], m);
753 scale_rows(mat, m, 1 + row);
754 isl_int_divexact(m, mat->row[0][0], tab->mat->row[r][0]);
755 isl_int_mul(mat->row[1 + row][0], m, tab->mat->row[r][1]);
756 for (j = 0; j < tab->n_param; ++j) {
757 int col;
758 if (tab->var[j].is_row)
759 continue;
760 col = tab->var[j].index;
761 isl_int_mul(mat->row[1 + row][1 + j], m,
762 tab->mat->row[r][off + col]);
764 for (j = 0; j < tab->n_div; ++j) {
765 int col;
766 if (tab->var[tab->n_var - tab->n_div+j].is_row)
767 continue;
768 col = tab->var[tab->n_var - tab->n_div+j].index;
769 isl_int_mul(mat->row[1 + row][1 + tab->n_param + j], m,
770 tab->mat->row[r][off + col]);
772 if (sol->max)
773 isl_seq_neg(mat->row[1 + row], mat->row[1 + row],
774 mat->n_col);
777 isl_int_clear(m);
779 sol_push_sol_mat(sol, bset, mat);
780 return;
781 error2:
782 isl_int_clear(m);
783 error:
784 isl_basic_set_free(bset);
785 isl_mat_free(mat);
786 sol->error = 1;
789 struct isl_sol_map {
790 struct isl_sol sol;
791 struct isl_map *map;
792 struct isl_set *empty;
795 static void sol_map_free(struct isl_sol *sol)
797 struct isl_sol_map *sol_map = (struct isl_sol_map *) sol;
798 isl_map_free(sol_map->map);
799 isl_set_free(sol_map->empty);
802 /* This function is called for parts of the context where there is
803 * no solution, with "bset" corresponding to the context tableau.
804 * Simply add the basic set to the set "empty".
806 static void sol_map_add_empty(struct isl_sol_map *sol,
807 struct isl_basic_set *bset)
809 if (!bset || !sol->empty)
810 goto error;
812 sol->empty = isl_set_grow(sol->empty, 1);
813 bset = isl_basic_set_simplify(bset);
814 bset = isl_basic_set_finalize(bset);
815 sol->empty = isl_set_add_basic_set(sol->empty, isl_basic_set_copy(bset));
816 if (!sol->empty)
817 goto error;
818 isl_basic_set_free(bset);
819 return;
820 error:
821 isl_basic_set_free(bset);
822 sol->sol.error = 1;
825 static void sol_map_add_empty_wrap(struct isl_sol *sol,
826 struct isl_basic_set *bset)
828 sol_map_add_empty((struct isl_sol_map *)sol, bset);
831 /* Given a basic set "dom" that represents the context and a tuple of
832 * affine expressions "ma" defined over this domain, construct a basic map
833 * that expresses this function on the domain.
835 static void sol_map_add(struct isl_sol_map *sol,
836 __isl_take isl_basic_set *dom, __isl_take isl_multi_aff *ma)
838 isl_basic_map *bmap;
840 if (sol->sol.error || !dom || !ma)
841 goto error;
843 bmap = isl_basic_map_from_multi_aff2(ma, sol->sol.rational);
844 bmap = isl_basic_map_intersect_domain(bmap, dom);
845 sol->map = isl_map_grow(sol->map, 1);
846 sol->map = isl_map_add_basic_map(sol->map, bmap);
847 if (!sol->map)
848 sol->sol.error = 1;
849 return;
850 error:
851 isl_basic_set_free(dom);
852 isl_multi_aff_free(ma);
853 sol->sol.error = 1;
856 static void sol_map_add_wrap(struct isl_sol *sol,
857 __isl_take isl_basic_set *dom, __isl_take isl_multi_aff *ma)
859 sol_map_add((struct isl_sol_map *)sol, dom, ma);
863 /* Store the "parametric constant" of row "row" of tableau "tab" in "line",
864 * i.e., the constant term and the coefficients of all variables that
865 * appear in the context tableau.
866 * Note that the coefficient of the big parameter M is NOT copied.
867 * The context tableau may not have a big parameter and even when it
868 * does, it is a different big parameter.
870 static void get_row_parameter_line(struct isl_tab *tab, int row, isl_int *line)
872 int i;
873 unsigned off = 2 + tab->M;
875 isl_int_set(line[0], tab->mat->row[row][1]);
876 for (i = 0; i < tab->n_param; ++i) {
877 if (tab->var[i].is_row)
878 isl_int_set_si(line[1 + i], 0);
879 else {
880 int col = tab->var[i].index;
881 isl_int_set(line[1 + i], tab->mat->row[row][off + col]);
884 for (i = 0; i < tab->n_div; ++i) {
885 if (tab->var[tab->n_var - tab->n_div + i].is_row)
886 isl_int_set_si(line[1 + tab->n_param + i], 0);
887 else {
888 int col = tab->var[tab->n_var - tab->n_div + i].index;
889 isl_int_set(line[1 + tab->n_param + i],
890 tab->mat->row[row][off + col]);
895 /* Check if rows "row1" and "row2" have identical "parametric constants",
896 * as explained above.
897 * In this case, we also insist that the coefficients of the big parameter
898 * be the same as the values of the constants will only be the same
899 * if these coefficients are also the same.
901 static int identical_parameter_line(struct isl_tab *tab, int row1, int row2)
903 int i;
904 unsigned off = 2 + tab->M;
906 if (isl_int_ne(tab->mat->row[row1][1], tab->mat->row[row2][1]))
907 return 0;
909 if (tab->M && isl_int_ne(tab->mat->row[row1][2],
910 tab->mat->row[row2][2]))
911 return 0;
913 for (i = 0; i < tab->n_param + tab->n_div; ++i) {
914 int pos = i < tab->n_param ? i :
915 tab->n_var - tab->n_div + i - tab->n_param;
916 int col;
918 if (tab->var[pos].is_row)
919 continue;
920 col = tab->var[pos].index;
921 if (isl_int_ne(tab->mat->row[row1][off + col],
922 tab->mat->row[row2][off + col]))
923 return 0;
925 return 1;
928 /* Return an inequality that expresses that the "parametric constant"
929 * should be non-negative.
930 * This function is only called when the coefficient of the big parameter
931 * is equal to zero.
933 static struct isl_vec *get_row_parameter_ineq(struct isl_tab *tab, int row)
935 struct isl_vec *ineq;
937 ineq = isl_vec_alloc(tab->mat->ctx, 1 + tab->n_param + tab->n_div);
938 if (!ineq)
939 return NULL;
941 get_row_parameter_line(tab, row, ineq->el);
942 if (ineq)
943 ineq = isl_vec_normalize(ineq);
945 return ineq;
948 /* Normalize a div expression of the form
950 * [(g*f(x) + c)/(g * m)]
952 * with c the constant term and f(x) the remaining coefficients, to
954 * [(f(x) + [c/g])/m]
956 static void normalize_div(__isl_keep isl_vec *div)
958 isl_ctx *ctx = isl_vec_get_ctx(div);
959 int len = div->size - 2;
961 isl_seq_gcd(div->el + 2, len, &ctx->normalize_gcd);
962 isl_int_gcd(ctx->normalize_gcd, ctx->normalize_gcd, div->el[0]);
964 if (isl_int_is_one(ctx->normalize_gcd))
965 return;
967 isl_int_divexact(div->el[0], div->el[0], ctx->normalize_gcd);
968 isl_int_fdiv_q(div->el[1], div->el[1], ctx->normalize_gcd);
969 isl_seq_scale_down(div->el + 2, div->el + 2, ctx->normalize_gcd, len);
972 /* Return an integer division for use in a parametric cut based
973 * on the given row.
974 * In particular, let the parametric constant of the row be
976 * \sum_i a_i y_i
978 * where y_0 = 1, but none of the y_i corresponds to the big parameter M.
979 * The div returned is equal to
981 * floor(\sum_i {-a_i} y_i) = floor((\sum_i (-a_i mod d) y_i)/d)
983 static struct isl_vec *get_row_parameter_div(struct isl_tab *tab, int row)
985 struct isl_vec *div;
987 div = isl_vec_alloc(tab->mat->ctx, 1 + 1 + tab->n_param + tab->n_div);
988 if (!div)
989 return NULL;
991 isl_int_set(div->el[0], tab->mat->row[row][0]);
992 get_row_parameter_line(tab, row, div->el + 1);
993 isl_seq_neg(div->el + 1, div->el + 1, div->size - 1);
994 normalize_div(div);
995 isl_seq_fdiv_r(div->el + 1, div->el + 1, div->el[0], div->size - 1);
997 return div;
1000 /* Return an integer division for use in transferring an integrality constraint
1001 * to the context.
1002 * In particular, let the parametric constant of the row be
1004 * \sum_i a_i y_i
1006 * where y_0 = 1, but none of the y_i corresponds to the big parameter M.
1007 * The the returned div is equal to
1009 * floor(\sum_i {a_i} y_i) = floor((\sum_i (a_i mod d) y_i)/d)
1011 static struct isl_vec *get_row_split_div(struct isl_tab *tab, int row)
1013 struct isl_vec *div;
1015 div = isl_vec_alloc(tab->mat->ctx, 1 + 1 + tab->n_param + tab->n_div);
1016 if (!div)
1017 return NULL;
1019 isl_int_set(div->el[0], tab->mat->row[row][0]);
1020 get_row_parameter_line(tab, row, div->el + 1);
1021 normalize_div(div);
1022 isl_seq_fdiv_r(div->el + 1, div->el + 1, div->el[0], div->size - 1);
1024 return div;
1027 /* Construct and return an inequality that expresses an upper bound
1028 * on the given div.
1029 * In particular, if the div is given by
1031 * d = floor(e/m)
1033 * then the inequality expresses
1035 * m d <= e
1037 static __isl_give isl_vec *ineq_for_div(__isl_keep isl_basic_set *bset,
1038 unsigned div)
1040 isl_size total;
1041 unsigned div_pos;
1042 struct isl_vec *ineq;
1044 total = isl_basic_set_dim(bset, isl_dim_all);
1045 if (total < 0)
1046 return NULL;
1048 div_pos = 1 + total - bset->n_div + div;
1050 ineq = isl_vec_alloc(bset->ctx, 1 + total);
1051 if (!ineq)
1052 return NULL;
1054 isl_seq_cpy(ineq->el, bset->div[div] + 1, 1 + total);
1055 isl_int_neg(ineq->el[div_pos], bset->div[div][0]);
1056 return ineq;
1059 /* Given a row in the tableau and a div that was created
1060 * using get_row_split_div and that has been constrained to equality, i.e.,
1062 * d = floor(\sum_i {a_i} y_i) = \sum_i {a_i} y_i
1064 * replace the expression "\sum_i {a_i} y_i" in the row by d,
1065 * i.e., we subtract "\sum_i {a_i} y_i" and add 1 d.
1066 * The coefficients of the non-parameters in the tableau have been
1067 * verified to be integral. We can therefore simply replace coefficient b
1068 * by floor(b). For the coefficients of the parameters we have
1069 * floor(a_i) = a_i - {a_i}, while for the other coefficients, we have
1070 * floor(b) = b.
1072 static struct isl_tab *set_row_cst_to_div(struct isl_tab *tab, int row, int div)
1074 isl_seq_fdiv_q(tab->mat->row[row] + 1, tab->mat->row[row] + 1,
1075 tab->mat->row[row][0], 1 + tab->M + tab->n_col);
1077 isl_int_set_si(tab->mat->row[row][0], 1);
1079 if (tab->var[tab->n_var - tab->n_div + div].is_row) {
1080 int drow = tab->var[tab->n_var - tab->n_div + div].index;
1082 isl_assert(tab->mat->ctx,
1083 isl_int_is_one(tab->mat->row[drow][0]), goto error);
1084 isl_seq_combine(tab->mat->row[row] + 1,
1085 tab->mat->ctx->one, tab->mat->row[row] + 1,
1086 tab->mat->ctx->one, tab->mat->row[drow] + 1,
1087 1 + tab->M + tab->n_col);
1088 } else {
1089 int dcol = tab->var[tab->n_var - tab->n_div + div].index;
1091 isl_int_add_ui(tab->mat->row[row][2 + tab->M + dcol],
1092 tab->mat->row[row][2 + tab->M + dcol], 1);
1095 return tab;
1096 error:
1097 isl_tab_free(tab);
1098 return NULL;
1101 /* Check if the (parametric) constant of the given row is obviously
1102 * negative, meaning that we don't need to consult the context tableau.
1103 * If there is a big parameter and its coefficient is non-zero,
1104 * then this coefficient determines the outcome.
1105 * Otherwise, we check whether the constant is negative and
1106 * all non-zero coefficients of parameters are negative and
1107 * belong to non-negative parameters.
1109 static int is_obviously_neg(struct isl_tab *tab, int row)
1111 int i;
1112 int col;
1113 unsigned off = 2 + tab->M;
1115 if (tab->M) {
1116 if (isl_int_is_pos(tab->mat->row[row][2]))
1117 return 0;
1118 if (isl_int_is_neg(tab->mat->row[row][2]))
1119 return 1;
1122 if (isl_int_is_nonneg(tab->mat->row[row][1]))
1123 return 0;
1124 for (i = 0; i < tab->n_param; ++i) {
1125 /* Eliminated parameter */
1126 if (tab->var[i].is_row)
1127 continue;
1128 col = tab->var[i].index;
1129 if (isl_int_is_zero(tab->mat->row[row][off + col]))
1130 continue;
1131 if (!tab->var[i].is_nonneg)
1132 return 0;
1133 if (isl_int_is_pos(tab->mat->row[row][off + col]))
1134 return 0;
1136 for (i = 0; i < tab->n_div; ++i) {
1137 if (tab->var[tab->n_var - tab->n_div + i].is_row)
1138 continue;
1139 col = tab->var[tab->n_var - tab->n_div + i].index;
1140 if (isl_int_is_zero(tab->mat->row[row][off + col]))
1141 continue;
1142 if (!tab->var[tab->n_var - tab->n_div + i].is_nonneg)
1143 return 0;
1144 if (isl_int_is_pos(tab->mat->row[row][off + col]))
1145 return 0;
1147 return 1;
1150 /* Check if the (parametric) constant of the given row is obviously
1151 * non-negative, meaning that we don't need to consult the context tableau.
1152 * If there is a big parameter and its coefficient is non-zero,
1153 * then this coefficient determines the outcome.
1154 * Otherwise, we check whether the constant is non-negative and
1155 * all non-zero coefficients of parameters are positive and
1156 * belong to non-negative parameters.
1158 static int is_obviously_nonneg(struct isl_tab *tab, int row)
1160 int i;
1161 int col;
1162 unsigned off = 2 + tab->M;
1164 if (tab->M) {
1165 if (isl_int_is_pos(tab->mat->row[row][2]))
1166 return 1;
1167 if (isl_int_is_neg(tab->mat->row[row][2]))
1168 return 0;
1171 if (isl_int_is_neg(tab->mat->row[row][1]))
1172 return 0;
1173 for (i = 0; i < tab->n_param; ++i) {
1174 /* Eliminated parameter */
1175 if (tab->var[i].is_row)
1176 continue;
1177 col = tab->var[i].index;
1178 if (isl_int_is_zero(tab->mat->row[row][off + col]))
1179 continue;
1180 if (!tab->var[i].is_nonneg)
1181 return 0;
1182 if (isl_int_is_neg(tab->mat->row[row][off + col]))
1183 return 0;
1185 for (i = 0; i < tab->n_div; ++i) {
1186 if (tab->var[tab->n_var - tab->n_div + i].is_row)
1187 continue;
1188 col = tab->var[tab->n_var - tab->n_div + i].index;
1189 if (isl_int_is_zero(tab->mat->row[row][off + col]))
1190 continue;
1191 if (!tab->var[tab->n_var - tab->n_div + i].is_nonneg)
1192 return 0;
1193 if (isl_int_is_neg(tab->mat->row[row][off + col]))
1194 return 0;
1196 return 1;
1199 /* Given a row r and two columns, return the column that would
1200 * lead to the lexicographically smallest increment in the sample
1201 * solution when leaving the basis in favor of the row.
1202 * Pivoting with column c will increment the sample value by a non-negative
1203 * constant times a_{V,c}/a_{r,c}, with a_{V,c} the elements of column c
1204 * corresponding to the non-parametric variables.
1205 * If variable v appears in a column c_v, then a_{v,c} = 1 iff c = c_v,
1206 * with all other entries in this virtual row equal to zero.
1207 * If variable v appears in a row, then a_{v,c} is the element in column c
1208 * of that row.
1210 * Let v be the first variable with a_{v,c1}/a_{r,c1} != a_{v,c2}/a_{r,c2}.
1211 * Then if a_{v,c1}/a_{r,c1} < a_{v,c2}/a_{r,c2}, i.e.,
1212 * a_{v,c2} a_{r,c1} - a_{v,c1} a_{r,c2} > 0, c1 results in the minimal
1213 * increment. Otherwise, it's c2.
1215 static int lexmin_col_pair(struct isl_tab *tab,
1216 int row, int col1, int col2, isl_int tmp)
1218 int i;
1219 isl_int *tr;
1221 tr = tab->mat->row[row] + 2 + tab->M;
1223 for (i = tab->n_param; i < tab->n_var - tab->n_div; ++i) {
1224 int s1, s2;
1225 isl_int *r;
1227 if (!tab->var[i].is_row) {
1228 if (tab->var[i].index == col1)
1229 return col2;
1230 if (tab->var[i].index == col2)
1231 return col1;
1232 continue;
1235 if (tab->var[i].index == row)
1236 continue;
1238 r = tab->mat->row[tab->var[i].index] + 2 + tab->M;
1239 s1 = isl_int_sgn(r[col1]);
1240 s2 = isl_int_sgn(r[col2]);
1241 if (s1 == 0 && s2 == 0)
1242 continue;
1243 if (s1 < s2)
1244 return col1;
1245 if (s2 < s1)
1246 return col2;
1248 isl_int_mul(tmp, r[col2], tr[col1]);
1249 isl_int_submul(tmp, r[col1], tr[col2]);
1250 if (isl_int_is_pos(tmp))
1251 return col1;
1252 if (isl_int_is_neg(tmp))
1253 return col2;
1255 return -1;
1258 /* Does the index into the tab->var or tab->con array "index"
1259 * correspond to a variable in the context tableau?
1260 * In particular, it needs to be an index into the tab->var array and
1261 * it needs to refer to either one of the first tab->n_param variables or
1262 * one of the last tab->n_div variables.
1264 static int is_parameter_var(struct isl_tab *tab, int index)
1266 if (index < 0)
1267 return 0;
1268 if (index < tab->n_param)
1269 return 1;
1270 if (index >= tab->n_var - tab->n_div)
1271 return 1;
1272 return 0;
1275 /* Does column "col" of "tab" refer to a variable in the context tableau?
1277 static int col_is_parameter_var(struct isl_tab *tab, int col)
1279 return is_parameter_var(tab, tab->col_var[col]);
1282 /* Does row "row" of "tab" refer to a variable in the context tableau?
1284 static int row_is_parameter_var(struct isl_tab *tab, int row)
1286 return is_parameter_var(tab, tab->row_var[row]);
1289 /* Given a row in the tableau, find and return the column that would
1290 * result in the lexicographically smallest, but positive, increment
1291 * in the sample point.
1292 * If there is no such column, then return tab->n_col.
1293 * If anything goes wrong, return -1.
1295 static int lexmin_pivot_col(struct isl_tab *tab, int row)
1297 int j;
1298 int col = tab->n_col;
1299 isl_int *tr;
1300 isl_int tmp;
1302 tr = tab->mat->row[row] + 2 + tab->M;
1304 isl_int_init(tmp);
1306 for (j = tab->n_dead; j < tab->n_col; ++j) {
1307 if (col_is_parameter_var(tab, j))
1308 continue;
1310 if (!isl_int_is_pos(tr[j]))
1311 continue;
1313 if (col == tab->n_col)
1314 col = j;
1315 else
1316 col = lexmin_col_pair(tab, row, col, j, tmp);
1317 isl_assert(tab->mat->ctx, col >= 0, goto error);
1320 isl_int_clear(tmp);
1321 return col;
1322 error:
1323 isl_int_clear(tmp);
1324 return -1;
1327 /* Return the first known violated constraint, i.e., a non-negative
1328 * constraint that currently has an either obviously negative value
1329 * or a previously determined to be negative value.
1331 * If any constraint has a negative coefficient for the big parameter,
1332 * if any, then we return one of these first.
1334 static int first_neg(struct isl_tab *tab)
1336 int row;
1338 if (tab->M)
1339 for (row = tab->n_redundant; row < tab->n_row; ++row) {
1340 if (!isl_tab_var_from_row(tab, row)->is_nonneg)
1341 continue;
1342 if (!isl_int_is_neg(tab->mat->row[row][2]))
1343 continue;
1344 if (tab->row_sign)
1345 tab->row_sign[row] = isl_tab_row_neg;
1346 return row;
1348 for (row = tab->n_redundant; row < tab->n_row; ++row) {
1349 if (!isl_tab_var_from_row(tab, row)->is_nonneg)
1350 continue;
1351 if (tab->row_sign) {
1352 if (tab->row_sign[row] == 0 &&
1353 is_obviously_neg(tab, row))
1354 tab->row_sign[row] = isl_tab_row_neg;
1355 if (tab->row_sign[row] != isl_tab_row_neg)
1356 continue;
1357 } else if (!is_obviously_neg(tab, row))
1358 continue;
1359 return row;
1361 return -1;
1364 /* Check whether the invariant that all columns are lexico-positive
1365 * is satisfied. This function is not called from the current code
1366 * but is useful during debugging.
1368 static void check_lexpos(struct isl_tab *tab) __attribute__ ((unused));
1369 static void check_lexpos(struct isl_tab *tab)
1371 unsigned off = 2 + tab->M;
1372 int col;
1373 int var;
1374 int row;
1376 for (col = tab->n_dead; col < tab->n_col; ++col) {
1377 if (col_is_parameter_var(tab, col))
1378 continue;
1379 for (var = tab->n_param; var < tab->n_var - tab->n_div; ++var) {
1380 if (!tab->var[var].is_row) {
1381 if (tab->var[var].index == col)
1382 break;
1383 else
1384 continue;
1386 row = tab->var[var].index;
1387 if (isl_int_is_zero(tab->mat->row[row][off + col]))
1388 continue;
1389 if (isl_int_is_pos(tab->mat->row[row][off + col]))
1390 break;
1391 fprintf(stderr, "lexneg column %d (row %d)\n",
1392 col, row);
1394 if (var >= tab->n_var - tab->n_div)
1395 fprintf(stderr, "zero column %d\n", col);
1399 /* Report to the caller that the given constraint is part of an encountered
1400 * conflict.
1402 static int report_conflicting_constraint(struct isl_tab *tab, int con)
1404 return tab->conflict(con, tab->conflict_user);
1407 /* Given a conflicting row in the tableau, report all constraints
1408 * involved in the row to the caller. That is, the row itself
1409 * (if it represents a constraint) and all constraint columns with
1410 * non-zero (and therefore negative) coefficients.
1412 static int report_conflict(struct isl_tab *tab, int row)
1414 int j;
1415 isl_int *tr;
1417 if (!tab->conflict)
1418 return 0;
1420 if (tab->row_var[row] < 0 &&
1421 report_conflicting_constraint(tab, ~tab->row_var[row]) < 0)
1422 return -1;
1424 tr = tab->mat->row[row] + 2 + tab->M;
1426 for (j = tab->n_dead; j < tab->n_col; ++j) {
1427 if (col_is_parameter_var(tab, j))
1428 continue;
1430 if (!isl_int_is_neg(tr[j]))
1431 continue;
1433 if (tab->col_var[j] < 0 &&
1434 report_conflicting_constraint(tab, ~tab->col_var[j]) < 0)
1435 return -1;
1438 return 0;
1441 /* Resolve all known or obviously violated constraints through pivoting.
1442 * In particular, as long as we can find any violated constraint, we
1443 * look for a pivoting column that would result in the lexicographically
1444 * smallest increment in the sample point. If there is no such column
1445 * then the tableau is infeasible.
1447 static int restore_lexmin(struct isl_tab *tab) WARN_UNUSED;
1448 static int restore_lexmin(struct isl_tab *tab)
1450 int row, col;
1452 if (!tab)
1453 return -1;
1454 if (tab->empty)
1455 return 0;
1456 while ((row = first_neg(tab)) != -1) {
1457 col = lexmin_pivot_col(tab, row);
1458 if (col >= tab->n_col) {
1459 if (report_conflict(tab, row) < 0)
1460 return -1;
1461 if (isl_tab_mark_empty(tab) < 0)
1462 return -1;
1463 return 0;
1465 if (col < 0)
1466 return -1;
1467 if (isl_tab_pivot(tab, row, col) < 0)
1468 return -1;
1470 return 0;
1473 /* Given a row that represents an equality, look for an appropriate
1474 * pivoting column.
1475 * In particular, if there are any non-zero coefficients among
1476 * the non-parameter variables, then we take the last of these
1477 * variables. Eliminating this variable in terms of the other
1478 * variables and/or parameters does not influence the property
1479 * that all column in the initial tableau are lexicographically
1480 * positive. The row corresponding to the eliminated variable
1481 * will only have non-zero entries below the diagonal of the
1482 * initial tableau. That is, we transform
1484 * I I
1485 * 1 into a
1486 * I I
1488 * If there is no such non-parameter variable, then we are dealing with
1489 * pure parameter equality and we pick any parameter with coefficient 1 or -1
1490 * for elimination. This will ensure that the eliminated parameter
1491 * always has an integer value whenever all the other parameters are integral.
1492 * If there is no such parameter then we return -1.
1494 static int last_var_col_or_int_par_col(struct isl_tab *tab, int row)
1496 unsigned off = 2 + tab->M;
1497 int i;
1499 for (i = tab->n_var - tab->n_div - 1; i >= 0 && i >= tab->n_param; --i) {
1500 int col;
1501 if (tab->var[i].is_row)
1502 continue;
1503 col = tab->var[i].index;
1504 if (col <= tab->n_dead)
1505 continue;
1506 if (!isl_int_is_zero(tab->mat->row[row][off + col]))
1507 return col;
1509 for (i = tab->n_dead; i < tab->n_col; ++i) {
1510 if (isl_int_is_one(tab->mat->row[row][off + i]))
1511 return i;
1512 if (isl_int_is_negone(tab->mat->row[row][off + i]))
1513 return i;
1515 return -1;
1518 /* Add an equality that is known to be valid to the tableau.
1519 * We first check if we can eliminate a variable or a parameter.
1520 * If not, we add the equality as two inequalities.
1521 * In this case, the equality was a pure parameter equality and there
1522 * is no need to resolve any constraint violations.
1524 * This function assumes that at least two more rows and at least
1525 * two more elements in the constraint array are available in the tableau.
1527 static struct isl_tab *add_lexmin_valid_eq(struct isl_tab *tab, isl_int *eq)
1529 int i;
1530 int r;
1532 if (!tab)
1533 return NULL;
1534 r = isl_tab_add_row(tab, eq);
1535 if (r < 0)
1536 goto error;
1538 r = tab->con[r].index;
1539 i = last_var_col_or_int_par_col(tab, r);
1540 if (i < 0) {
1541 tab->con[r].is_nonneg = 1;
1542 if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r]) < 0)
1543 goto error;
1544 isl_seq_neg(eq, eq, 1 + tab->n_var);
1545 r = isl_tab_add_row(tab, eq);
1546 if (r < 0)
1547 goto error;
1548 tab->con[r].is_nonneg = 1;
1549 if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r]) < 0)
1550 goto error;
1551 } else {
1552 if (isl_tab_pivot(tab, r, i) < 0)
1553 goto error;
1554 if (isl_tab_kill_col(tab, i) < 0)
1555 goto error;
1556 tab->n_eq++;
1559 return tab;
1560 error:
1561 isl_tab_free(tab);
1562 return NULL;
1565 /* Check if the given row is a pure constant.
1567 static int is_constant(struct isl_tab *tab, int row)
1569 unsigned off = 2 + tab->M;
1571 return !isl_seq_any_non_zero(tab->mat->row[row] + off + tab->n_dead,
1572 tab->n_col - tab->n_dead);
1575 /* Is the given row a parametric constant?
1576 * That is, does it only involve variables that also appear in the context?
1578 static int is_parametric_constant(struct isl_tab *tab, int row)
1580 unsigned off = 2 + tab->M;
1581 int col;
1583 for (col = tab->n_dead; col < tab->n_col; ++col) {
1584 if (col_is_parameter_var(tab, col))
1585 continue;
1586 if (isl_int_is_zero(tab->mat->row[row][off + col]))
1587 continue;
1588 return 0;
1591 return 1;
1594 /* Add an equality that may or may not be valid to the tableau.
1595 * If the resulting row is a pure constant, then it must be zero.
1596 * Otherwise, the resulting tableau is empty.
1598 * If the row is not a pure constant, then we add two inequalities,
1599 * each time checking that they can be satisfied.
1600 * In the end we try to use one of the two constraints to eliminate
1601 * a column.
1603 * This function assumes that at least two more rows and at least
1604 * two more elements in the constraint array are available in the tableau.
1606 static int add_lexmin_eq(struct isl_tab *tab, isl_int *eq) WARN_UNUSED;
1607 static int add_lexmin_eq(struct isl_tab *tab, isl_int *eq)
1609 int r1, r2;
1610 int row;
1611 struct isl_tab_undo *snap;
1613 if (!tab)
1614 return -1;
1615 snap = isl_tab_snap(tab);
1616 r1 = isl_tab_add_row(tab, eq);
1617 if (r1 < 0)
1618 return -1;
1619 tab->con[r1].is_nonneg = 1;
1620 if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r1]) < 0)
1621 return -1;
1623 row = tab->con[r1].index;
1624 if (is_constant(tab, row)) {
1625 if (!isl_int_is_zero(tab->mat->row[row][1]) ||
1626 (tab->M && !isl_int_is_zero(tab->mat->row[row][2]))) {
1627 if (isl_tab_mark_empty(tab) < 0)
1628 return -1;
1629 return 0;
1631 if (isl_tab_rollback(tab, snap) < 0)
1632 return -1;
1633 return 0;
1636 if (restore_lexmin(tab) < 0)
1637 return -1;
1638 if (tab->empty)
1639 return 0;
1641 isl_seq_neg(eq, eq, 1 + tab->n_var);
1643 r2 = isl_tab_add_row(tab, eq);
1644 if (r2 < 0)
1645 return -1;
1646 tab->con[r2].is_nonneg = 1;
1647 if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r2]) < 0)
1648 return -1;
1650 if (restore_lexmin(tab) < 0)
1651 return -1;
1652 if (tab->empty)
1653 return 0;
1655 if (!tab->con[r1].is_row) {
1656 if (isl_tab_kill_col(tab, tab->con[r1].index) < 0)
1657 return -1;
1658 } else if (!tab->con[r2].is_row) {
1659 if (isl_tab_kill_col(tab, tab->con[r2].index) < 0)
1660 return -1;
1663 if (tab->bmap) {
1664 tab->bmap = isl_basic_map_add_ineq(tab->bmap, eq);
1665 if (isl_tab_push(tab, isl_tab_undo_bmap_ineq) < 0)
1666 return -1;
1667 isl_seq_neg(eq, eq, 1 + tab->n_var);
1668 tab->bmap = isl_basic_map_add_ineq(tab->bmap, eq);
1669 isl_seq_neg(eq, eq, 1 + tab->n_var);
1670 if (isl_tab_push(tab, isl_tab_undo_bmap_ineq) < 0)
1671 return -1;
1672 if (!tab->bmap)
1673 return -1;
1676 return 0;
1679 /* Add an inequality to the tableau, resolving violations using
1680 * restore_lexmin.
1682 * This function assumes that at least one more row and at least
1683 * one more element in the constraint array are available in the tableau.
1685 static struct isl_tab *add_lexmin_ineq(struct isl_tab *tab, isl_int *ineq)
1687 int r;
1689 if (!tab)
1690 return NULL;
1691 if (tab->bmap) {
1692 tab->bmap = isl_basic_map_add_ineq(tab->bmap, ineq);
1693 if (isl_tab_push(tab, isl_tab_undo_bmap_ineq) < 0)
1694 goto error;
1695 if (!tab->bmap)
1696 goto error;
1698 r = isl_tab_add_row(tab, ineq);
1699 if (r < 0)
1700 goto error;
1701 tab->con[r].is_nonneg = 1;
1702 if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r]) < 0)
1703 goto error;
1704 if (isl_tab_row_is_redundant(tab, tab->con[r].index)) {
1705 if (isl_tab_mark_redundant(tab, tab->con[r].index) < 0)
1706 goto error;
1707 return tab;
1710 if (restore_lexmin(tab) < 0)
1711 goto error;
1712 if (!tab->empty && tab->con[r].is_row &&
1713 isl_tab_row_is_redundant(tab, tab->con[r].index))
1714 if (isl_tab_mark_redundant(tab, tab->con[r].index) < 0)
1715 goto error;
1716 return tab;
1717 error:
1718 isl_tab_free(tab);
1719 return NULL;
1722 /* Check if the coefficients of the parameters are all integral.
1724 static int integer_parameter(struct isl_tab *tab, int row)
1726 int i;
1727 int col;
1728 unsigned off = 2 + tab->M;
1730 for (i = 0; i < tab->n_param; ++i) {
1731 /* Eliminated parameter */
1732 if (tab->var[i].is_row)
1733 continue;
1734 col = tab->var[i].index;
1735 if (!isl_int_is_divisible_by(tab->mat->row[row][off + col],
1736 tab->mat->row[row][0]))
1737 return 0;
1739 for (i = 0; i < tab->n_div; ++i) {
1740 if (tab->var[tab->n_var - tab->n_div + i].is_row)
1741 continue;
1742 col = tab->var[tab->n_var - tab->n_div + i].index;
1743 if (!isl_int_is_divisible_by(tab->mat->row[row][off + col],
1744 tab->mat->row[row][0]))
1745 return 0;
1747 return 1;
1750 /* Check if the coefficients of the non-parameter variables are all integral.
1752 static int integer_variable(struct isl_tab *tab, int row)
1754 int i;
1755 unsigned off = 2 + tab->M;
1757 for (i = tab->n_dead; i < tab->n_col; ++i) {
1758 if (col_is_parameter_var(tab, i))
1759 continue;
1760 if (!isl_int_is_divisible_by(tab->mat->row[row][off + i],
1761 tab->mat->row[row][0]))
1762 return 0;
1764 return 1;
1767 /* Check if the constant term is integral.
1769 static int integer_constant(struct isl_tab *tab, int row)
1771 return isl_int_is_divisible_by(tab->mat->row[row][1],
1772 tab->mat->row[row][0]);
1775 #define I_CST 1 << 0
1776 #define I_PAR 1 << 1
1777 #define I_VAR 1 << 2
1779 /* Check for next (non-parameter) variable after "var" (first if var == -1)
1780 * that is non-integer and therefore requires a cut and return
1781 * the index of the variable.
1782 * For parametric tableaus, there are three parts in a row,
1783 * the constant, the coefficients of the parameters and the rest.
1784 * For each part, we check whether the coefficients in that part
1785 * are all integral and if so, set the corresponding flag in *f.
1786 * If the constant and the parameter part are integral, then the
1787 * current sample value is integral and no cut is required
1788 * (irrespective of whether the variable part is integral).
1790 static int next_non_integer_var(struct isl_tab *tab, int var, int *f)
1792 var = var < 0 ? tab->n_param : var + 1;
1794 for (; var < tab->n_var - tab->n_div; ++var) {
1795 int flags = 0;
1796 int row;
1797 if (!tab->var[var].is_row)
1798 continue;
1799 row = tab->var[var].index;
1800 if (integer_constant(tab, row))
1801 ISL_FL_SET(flags, I_CST);
1802 if (integer_parameter(tab, row))
1803 ISL_FL_SET(flags, I_PAR);
1804 if (ISL_FL_ISSET(flags, I_CST) && ISL_FL_ISSET(flags, I_PAR))
1805 continue;
1806 if (integer_variable(tab, row))
1807 ISL_FL_SET(flags, I_VAR);
1808 *f = flags;
1809 return var;
1811 return -1;
1814 /* Check for first (non-parameter) variable that is non-integer and
1815 * therefore requires a cut and return the corresponding row.
1816 * For parametric tableaus, there are three parts in a row,
1817 * the constant, the coefficients of the parameters and the rest.
1818 * For each part, we check whether the coefficients in that part
1819 * are all integral and if so, set the corresponding flag in *f.
1820 * If the constant and the parameter part are integral, then the
1821 * current sample value is integral and no cut is required
1822 * (irrespective of whether the variable part is integral).
1824 static int first_non_integer_row(struct isl_tab *tab, int *f)
1826 int var = next_non_integer_var(tab, -1, f);
1828 return var < 0 ? -1 : tab->var[var].index;
1831 /* Add a (non-parametric) cut to cut away the non-integral sample
1832 * value of the given row.
1834 * If the row is given by
1836 * m r = f + \sum_i a_i y_i
1838 * then the cut is
1840 * c = - {-f/m} + \sum_i {a_i/m} y_i >= 0
1842 * The big parameter, if any, is ignored, since it is assumed to be big
1843 * enough to be divisible by any integer.
1844 * If the tableau is actually a parametric tableau, then this function
1845 * is only called when all coefficients of the parameters are integral.
1846 * The cut therefore has zero coefficients for the parameters.
1848 * The current value is known to be negative, so row_sign, if it
1849 * exists, is set accordingly.
1851 * Return the row of the cut or -1.
1853 static int add_cut(struct isl_tab *tab, int row)
1855 int i;
1856 int r;
1857 isl_int *r_row;
1858 unsigned off = 2 + tab->M;
1860 if (isl_tab_extend_cons(tab, 1) < 0)
1861 return -1;
1862 r = isl_tab_allocate_con(tab);
1863 if (r < 0)
1864 return -1;
1866 r_row = tab->mat->row[tab->con[r].index];
1867 isl_int_set(r_row[0], tab->mat->row[row][0]);
1868 isl_int_neg(r_row[1], tab->mat->row[row][1]);
1869 isl_int_fdiv_r(r_row[1], r_row[1], tab->mat->row[row][0]);
1870 isl_int_neg(r_row[1], r_row[1]);
1871 if (tab->M)
1872 isl_int_set_si(r_row[2], 0);
1873 for (i = 0; i < tab->n_col; ++i)
1874 isl_int_fdiv_r(r_row[off + i],
1875 tab->mat->row[row][off + i], tab->mat->row[row][0]);
1877 tab->con[r].is_nonneg = 1;
1878 if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r]) < 0)
1879 return -1;
1880 if (tab->row_sign)
1881 tab->row_sign[tab->con[r].index] = isl_tab_row_neg;
1883 return tab->con[r].index;
1886 #define CUT_ALL 1
1887 #define CUT_ONE 0
1889 /* Given a non-parametric tableau, add cuts until an integer
1890 * sample point is obtained or until the tableau is determined
1891 * to be integer infeasible.
1892 * As long as there is any non-integer value in the sample point,
1893 * we add appropriate cuts, if possible, for each of these
1894 * non-integer values and then resolve the violated
1895 * cut constraints using restore_lexmin.
1896 * If one of the corresponding rows is equal to an integral
1897 * combination of variables/constraints plus a non-integral constant,
1898 * then there is no way to obtain an integer point and we return
1899 * a tableau that is marked empty.
1900 * The parameter cutting_strategy controls the strategy used when adding cuts
1901 * to remove non-integer points. CUT_ALL adds all possible cuts
1902 * before continuing the search. CUT_ONE adds only one cut at a time.
1904 static struct isl_tab *cut_to_integer_lexmin(struct isl_tab *tab,
1905 int cutting_strategy)
1907 int var;
1908 int row;
1909 int flags;
1911 if (!tab)
1912 return NULL;
1913 if (tab->empty)
1914 return tab;
1916 while ((var = next_non_integer_var(tab, -1, &flags)) != -1) {
1917 do {
1918 if (ISL_FL_ISSET(flags, I_VAR)) {
1919 if (isl_tab_mark_empty(tab) < 0)
1920 goto error;
1921 return tab;
1923 row = tab->var[var].index;
1924 row = add_cut(tab, row);
1925 if (row < 0)
1926 goto error;
1927 if (cutting_strategy == CUT_ONE)
1928 break;
1929 } while ((var = next_non_integer_var(tab, var, &flags)) != -1);
1930 if (restore_lexmin(tab) < 0)
1931 goto error;
1932 if (tab->empty)
1933 break;
1935 return tab;
1936 error:
1937 isl_tab_free(tab);
1938 return NULL;
1941 /* Check whether all the currently active samples also satisfy the inequality
1942 * "ineq" (treated as an equality if eq is set).
1943 * Remove those samples that do not.
1945 static struct isl_tab *check_samples(struct isl_tab *tab, isl_int *ineq, int eq)
1947 int i;
1948 isl_int v;
1950 if (!tab)
1951 return NULL;
1953 isl_assert(tab->mat->ctx, tab->bmap, goto error);
1954 isl_assert(tab->mat->ctx, tab->samples, goto error);
1955 isl_assert(tab->mat->ctx, tab->samples->n_col == 1 + tab->n_var, goto error);
1957 isl_int_init(v);
1958 for (i = tab->n_outside; i < tab->n_sample; ++i) {
1959 int sgn;
1960 isl_seq_inner_product(ineq, tab->samples->row[i],
1961 1 + tab->n_var, &v);
1962 sgn = isl_int_sgn(v);
1963 if (eq ? (sgn == 0) : (sgn >= 0))
1964 continue;
1965 tab = isl_tab_drop_sample(tab, i);
1966 if (!tab)
1967 break;
1969 isl_int_clear(v);
1971 return tab;
1972 error:
1973 isl_tab_free(tab);
1974 return NULL;
1977 /* Check whether the sample value of the tableau is finite,
1978 * i.e., either the tableau does not use a big parameter, or
1979 * all values of the variables are equal to the big parameter plus
1980 * some constant. This constant is the actual sample value.
1982 static int sample_is_finite(struct isl_tab *tab)
1984 int i;
1986 if (!tab->M)
1987 return 1;
1989 for (i = 0; i < tab->n_var; ++i) {
1990 int row;
1991 if (!tab->var[i].is_row)
1992 return 0;
1993 row = tab->var[i].index;
1994 if (isl_int_ne(tab->mat->row[row][0], tab->mat->row[row][2]))
1995 return 0;
1997 return 1;
2000 /* Check if the context tableau of sol has any integer points.
2001 * Leave tab in empty state if no integer point can be found.
2002 * If an integer point can be found and if moreover it is finite,
2003 * then it is added to the list of sample values.
2005 * This function is only called when none of the currently active sample
2006 * values satisfies the most recently added constraint.
2008 static struct isl_tab *check_integer_feasible(struct isl_tab *tab)
2010 struct isl_tab_undo *snap;
2012 if (!tab)
2013 return NULL;
2015 snap = isl_tab_snap(tab);
2016 if (isl_tab_push_basis(tab) < 0)
2017 goto error;
2019 tab = cut_to_integer_lexmin(tab, CUT_ALL);
2020 if (!tab)
2021 goto error;
2023 if (!tab->empty && sample_is_finite(tab)) {
2024 struct isl_vec *sample;
2026 sample = isl_tab_get_sample_value(tab);
2028 if (isl_tab_add_sample(tab, sample) < 0)
2029 goto error;
2032 if (!tab->empty && isl_tab_rollback(tab, snap) < 0)
2033 goto error;
2035 return tab;
2036 error:
2037 isl_tab_free(tab);
2038 return NULL;
2041 /* Check if any of the currently active sample values satisfies
2042 * the inequality "ineq" (an equality if eq is set).
2044 static int tab_has_valid_sample(struct isl_tab *tab, isl_int *ineq, int eq)
2046 int i;
2047 isl_int v;
2049 if (!tab)
2050 return -1;
2052 isl_assert(tab->mat->ctx, tab->bmap, return -1);
2053 isl_assert(tab->mat->ctx, tab->samples, return -1);
2054 isl_assert(tab->mat->ctx, tab->samples->n_col == 1 + tab->n_var, return -1);
2056 isl_int_init(v);
2057 for (i = tab->n_outside; i < tab->n_sample; ++i) {
2058 int sgn;
2059 isl_seq_inner_product(ineq, tab->samples->row[i],
2060 1 + tab->n_var, &v);
2061 sgn = isl_int_sgn(v);
2062 if (eq ? (sgn == 0) : (sgn >= 0))
2063 break;
2065 isl_int_clear(v);
2067 return i < tab->n_sample;
2070 /* Insert a div specified by "div" to the tableau "tab" at position "pos" and
2071 * return isl_bool_true if the div is obviously non-negative.
2073 static isl_bool context_tab_insert_div(struct isl_tab *tab, int pos,
2074 __isl_keep isl_vec *div,
2075 isl_stat (*add_ineq)(void *user, isl_int *), void *user)
2077 int i;
2078 int r;
2079 struct isl_mat *samples;
2080 int nonneg;
2082 r = isl_tab_insert_div(tab, pos, div, add_ineq, user);
2083 if (r < 0)
2084 return isl_bool_error;
2085 nonneg = tab->var[r].is_nonneg;
2086 tab->var[r].frozen = 1;
2088 samples = isl_mat_extend(tab->samples,
2089 tab->n_sample, 1 + tab->n_var);
2090 tab->samples = samples;
2091 if (!samples)
2092 return isl_bool_error;
2093 for (i = tab->n_outside; i < samples->n_row; ++i) {
2094 isl_seq_inner_product(div->el + 1, samples->row[i],
2095 div->size - 1, &samples->row[i][samples->n_col - 1]);
2096 isl_int_fdiv_q(samples->row[i][samples->n_col - 1],
2097 samples->row[i][samples->n_col - 1], div->el[0]);
2099 tab->samples = isl_mat_move_cols(tab->samples, 1 + pos,
2100 1 + tab->n_var - 1, 1);
2101 if (!tab->samples)
2102 return isl_bool_error;
2104 return isl_bool_ok(nonneg);
2107 /* Add a div specified by "div" to both the main tableau and
2108 * the context tableau. In case of the main tableau, we only
2109 * need to add an extra div. In the context tableau, we also
2110 * need to express the meaning of the div.
2111 * Return the index of the div or -1 if anything went wrong.
2113 * The new integer division is added before any unknown integer
2114 * divisions in the context to ensure that it does not get
2115 * equated to some linear combination involving unknown integer
2116 * divisions.
2118 static int add_div(struct isl_tab *tab, struct isl_context *context,
2119 __isl_keep isl_vec *div)
2121 int r;
2122 int pos;
2123 isl_bool nonneg;
2124 struct isl_tab *context_tab = context->op->peek_tab(context);
2126 if (!tab || !context_tab)
2127 goto error;
2129 pos = context_tab->n_var - context->n_unknown;
2130 if ((nonneg = context->op->insert_div(context, pos, div)) < 0)
2131 goto error;
2133 if (!context->op->is_ok(context))
2134 goto error;
2136 pos = tab->n_var - context->n_unknown;
2137 if (isl_tab_extend_vars(tab, 1) < 0)
2138 goto error;
2139 r = isl_tab_insert_var(tab, pos);
2140 if (r < 0)
2141 goto error;
2142 if (nonneg)
2143 tab->var[r].is_nonneg = 1;
2144 tab->var[r].frozen = 1;
2145 tab->n_div++;
2147 return tab->n_div - 1 - context->n_unknown;
2148 error:
2149 context->op->invalidate(context);
2150 return -1;
2153 /* Return the position of the integer division that is equal to div/denom
2154 * if there is one. Otherwise, return a position beyond the integer divisions.
2156 static int find_div(struct isl_tab *tab, isl_int *div, isl_int denom)
2158 int i;
2159 isl_size total = isl_basic_map_dim(tab->bmap, isl_dim_all);
2160 isl_size n_div;
2162 n_div = isl_basic_map_dim(tab->bmap, isl_dim_div);
2163 if (total < 0 || n_div < 0)
2164 return -1;
2165 for (i = 0; i < n_div; ++i) {
2166 if (isl_int_ne(tab->bmap->div[i][0], denom))
2167 continue;
2168 if (!isl_seq_eq(tab->bmap->div[i] + 1, div, 1 + total))
2169 continue;
2170 return i;
2172 return n_div;
2175 /* Return the index of a div that corresponds to "div".
2176 * We first check if we already have such a div and if not, we create one.
2178 static int get_div(struct isl_tab *tab, struct isl_context *context,
2179 struct isl_vec *div)
2181 int d;
2182 struct isl_tab *context_tab = context->op->peek_tab(context);
2183 unsigned n_div;
2185 if (!context_tab)
2186 return -1;
2188 n_div = isl_basic_map_dim(context_tab->bmap, isl_dim_div);
2189 d = find_div(context_tab, div->el + 1, div->el[0]);
2190 if (d < 0)
2191 return -1;
2192 if (d < n_div)
2193 return d;
2195 return add_div(tab, context, div);
2198 /* Add a parametric cut to cut away the non-integral sample value
2199 * of the given row.
2200 * Let a_i be the coefficients of the constant term and the parameters
2201 * and let b_i be the coefficients of the variables or constraints
2202 * in basis of the tableau.
2203 * Let q be the div q = floor(\sum_i {-a_i} y_i).
2205 * The cut is expressed as
2207 * c = \sum_i -{-a_i} y_i + \sum_i {b_i} x_i + q >= 0
2209 * If q did not already exist in the context tableau, then it is added first.
2210 * If q is in a column of the main tableau then the "+ q" can be accomplished
2211 * by setting the corresponding entry to the denominator of the constraint.
2212 * If q happens to be in a row of the main tableau, then the corresponding
2213 * row needs to be added instead (taking care of the denominators).
2214 * Note that this is very unlikely, but perhaps not entirely impossible.
2216 * The current value of the cut is known to be negative (or at least
2217 * non-positive), so row_sign is set accordingly.
2219 * Return the row of the cut or -1.
2221 static int add_parametric_cut(struct isl_tab *tab, int row,
2222 struct isl_context *context)
2224 struct isl_vec *div;
2225 int d;
2226 int i;
2227 int r;
2228 isl_int *r_row;
2229 int col;
2230 int n;
2231 unsigned off = 2 + tab->M;
2233 if (!context)
2234 return -1;
2236 div = get_row_parameter_div(tab, row);
2237 if (!div)
2238 return -1;
2240 n = tab->n_div - context->n_unknown;
2241 d = context->op->get_div(context, tab, div);
2242 isl_vec_free(div);
2243 if (d < 0)
2244 return -1;
2246 if (isl_tab_extend_cons(tab, 1) < 0)
2247 return -1;
2248 r = isl_tab_allocate_con(tab);
2249 if (r < 0)
2250 return -1;
2252 r_row = tab->mat->row[tab->con[r].index];
2253 isl_int_set(r_row[0], tab->mat->row[row][0]);
2254 isl_int_neg(r_row[1], tab->mat->row[row][1]);
2255 isl_int_fdiv_r(r_row[1], r_row[1], tab->mat->row[row][0]);
2256 isl_int_neg(r_row[1], r_row[1]);
2257 if (tab->M)
2258 isl_int_set_si(r_row[2], 0);
2259 for (i = 0; i < tab->n_param; ++i) {
2260 if (tab->var[i].is_row)
2261 continue;
2262 col = tab->var[i].index;
2263 isl_int_neg(r_row[off + col], tab->mat->row[row][off + col]);
2264 isl_int_fdiv_r(r_row[off + col], r_row[off + col],
2265 tab->mat->row[row][0]);
2266 isl_int_neg(r_row[off + col], r_row[off + col]);
2268 for (i = 0; i < tab->n_div; ++i) {
2269 if (tab->var[tab->n_var - tab->n_div + i].is_row)
2270 continue;
2271 col = tab->var[tab->n_var - tab->n_div + i].index;
2272 isl_int_neg(r_row[off + col], tab->mat->row[row][off + col]);
2273 isl_int_fdiv_r(r_row[off + col], r_row[off + col],
2274 tab->mat->row[row][0]);
2275 isl_int_neg(r_row[off + col], r_row[off + col]);
2277 for (i = 0; i < tab->n_col; ++i) {
2278 if (tab->col_var[i] >= 0 &&
2279 (tab->col_var[i] < tab->n_param ||
2280 tab->col_var[i] >= tab->n_var - tab->n_div))
2281 continue;
2282 isl_int_fdiv_r(r_row[off + i],
2283 tab->mat->row[row][off + i], tab->mat->row[row][0]);
2285 if (tab->var[tab->n_var - tab->n_div + d].is_row) {
2286 isl_int gcd;
2287 int d_row = tab->var[tab->n_var - tab->n_div + d].index;
2288 isl_int_init(gcd);
2289 isl_int_gcd(gcd, tab->mat->row[d_row][0], r_row[0]);
2290 isl_int_divexact(r_row[0], r_row[0], gcd);
2291 isl_int_divexact(gcd, tab->mat->row[d_row][0], gcd);
2292 isl_seq_combine(r_row + 1, gcd, r_row + 1,
2293 r_row[0], tab->mat->row[d_row] + 1,
2294 off - 1 + tab->n_col);
2295 isl_int_mul(r_row[0], r_row[0], tab->mat->row[d_row][0]);
2296 isl_int_clear(gcd);
2297 } else {
2298 col = tab->var[tab->n_var - tab->n_div + d].index;
2299 isl_int_set(r_row[off + col], tab->mat->row[row][0]);
2302 tab->con[r].is_nonneg = 1;
2303 if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r]) < 0)
2304 return -1;
2305 if (tab->row_sign)
2306 tab->row_sign[tab->con[r].index] = isl_tab_row_neg;
2308 row = tab->con[r].index;
2310 if (d >= n && context->op->detect_equalities(context, tab) < 0)
2311 return -1;
2313 return row;
2316 /* Construct a tableau for bmap that can be used for computing
2317 * the lexicographic minimum (or maximum) of bmap.
2318 * If not NULL, then dom is the domain where the minimum
2319 * should be computed. In this case, we set up a parametric
2320 * tableau with row signs (initialized to "unknown").
2321 * If M is set, then the tableau will use a big parameter.
2322 * If max is set, then a maximum should be computed instead of a minimum.
2323 * This means that for each variable x, the tableau will contain the variable
2324 * x' = M - x, rather than x' = M + x. This in turn means that the coefficient
2325 * of the variables in all constraints are negated prior to adding them
2326 * to the tableau.
2328 static __isl_give struct isl_tab *tab_for_lexmin(__isl_keep isl_basic_map *bmap,
2329 __isl_keep isl_basic_set *dom, unsigned M, int max)
2331 int i;
2332 struct isl_tab *tab;
2333 unsigned n_var;
2334 unsigned o_var;
2335 isl_size total;
2337 total = isl_basic_map_dim(bmap, isl_dim_all);
2338 if (total < 0)
2339 return NULL;
2340 tab = isl_tab_alloc(bmap->ctx, 2 * bmap->n_eq + bmap->n_ineq + 1,
2341 total, M);
2342 if (!tab)
2343 return NULL;
2345 tab->rational = ISL_F_ISSET(bmap, ISL_BASIC_MAP_RATIONAL);
2346 if (dom) {
2347 isl_size dom_total;
2348 dom_total = isl_basic_set_dim(dom, isl_dim_all);
2349 if (dom_total < 0)
2350 goto error;
2351 tab->n_param = dom_total - dom->n_div;
2352 tab->n_div = dom->n_div;
2353 tab->row_sign = isl_calloc_array(bmap->ctx,
2354 enum isl_tab_row_sign, tab->mat->n_row);
2355 if (tab->mat->n_row && !tab->row_sign)
2356 goto error;
2358 if (ISL_F_ISSET(bmap, ISL_BASIC_MAP_EMPTY)) {
2359 if (isl_tab_mark_empty(tab) < 0)
2360 goto error;
2361 return tab;
2364 for (i = tab->n_param; i < tab->n_var - tab->n_div; ++i) {
2365 tab->var[i].is_nonneg = 1;
2366 tab->var[i].frozen = 1;
2368 o_var = 1 + tab->n_param;
2369 n_var = tab->n_var - tab->n_param - tab->n_div;
2370 for (i = 0; i < bmap->n_eq; ++i) {
2371 if (max)
2372 isl_seq_neg(bmap->eq[i] + o_var,
2373 bmap->eq[i] + o_var, n_var);
2374 tab = add_lexmin_valid_eq(tab, bmap->eq[i]);
2375 if (max)
2376 isl_seq_neg(bmap->eq[i] + o_var,
2377 bmap->eq[i] + o_var, n_var);
2378 if (!tab || tab->empty)
2379 return tab;
2381 if (bmap->n_eq && restore_lexmin(tab) < 0)
2382 goto error;
2383 for (i = 0; i < bmap->n_ineq; ++i) {
2384 if (max)
2385 isl_seq_neg(bmap->ineq[i] + o_var,
2386 bmap->ineq[i] + o_var, n_var);
2387 tab = add_lexmin_ineq(tab, bmap->ineq[i]);
2388 if (max)
2389 isl_seq_neg(bmap->ineq[i] + o_var,
2390 bmap->ineq[i] + o_var, n_var);
2391 if (!tab || tab->empty)
2392 return tab;
2394 return tab;
2395 error:
2396 isl_tab_free(tab);
2397 return NULL;
2400 /* Given a main tableau where more than one row requires a split,
2401 * determine and return the "best" row to split on.
2403 * If any of the rows requiring a split only involves
2404 * variables that also appear in the context tableau,
2405 * then the negative part is guaranteed not to have a solution.
2406 * It is therefore best to split on any of these rows first.
2408 * Otherwise,
2409 * given two rows in the main tableau, if the inequality corresponding
2410 * to the first row is redundant with respect to that of the second row
2411 * in the current tableau, then it is better to split on the second row,
2412 * since in the positive part, both rows will be positive.
2413 * (In the negative part a pivot will have to be performed and just about
2414 * anything can happen to the sign of the other row.)
2416 * As a simple heuristic, we therefore select the row that makes the most
2417 * of the other rows redundant.
2419 * Perhaps it would also be useful to look at the number of constraints
2420 * that conflict with any given constraint.
2422 * best is the best row so far (-1 when we have not found any row yet).
2423 * best_r is the number of other rows made redundant by row best.
2424 * When best is still -1, bset_r is meaningless, but it is initialized
2425 * to some arbitrary value (0) anyway. Without this redundant initialization
2426 * valgrind may warn about uninitialized memory accesses when isl
2427 * is compiled with some versions of gcc.
2429 static int best_split(struct isl_tab *tab, struct isl_tab *context_tab)
2431 struct isl_tab_undo *snap;
2432 int split;
2433 int row;
2434 int best = -1;
2435 int best_r = 0;
2437 if (isl_tab_extend_cons(context_tab, 2) < 0)
2438 return -1;
2440 snap = isl_tab_snap(context_tab);
2442 for (split = tab->n_redundant; split < tab->n_row; ++split) {
2443 struct isl_tab_undo *snap2;
2444 struct isl_vec *ineq = NULL;
2445 int r = 0;
2446 int ok;
2448 if (!isl_tab_var_from_row(tab, split)->is_nonneg)
2449 continue;
2450 if (tab->row_sign[split] != isl_tab_row_any)
2451 continue;
2453 if (is_parametric_constant(tab, split))
2454 return split;
2456 ineq = get_row_parameter_ineq(tab, split);
2457 if (!ineq)
2458 return -1;
2459 ok = isl_tab_add_ineq(context_tab, ineq->el) >= 0;
2460 isl_vec_free(ineq);
2461 if (!ok)
2462 return -1;
2464 snap2 = isl_tab_snap(context_tab);
2466 for (row = tab->n_redundant; row < tab->n_row; ++row) {
2467 struct isl_tab_var *var;
2469 if (row == split)
2470 continue;
2471 if (!isl_tab_var_from_row(tab, row)->is_nonneg)
2472 continue;
2473 if (tab->row_sign[row] != isl_tab_row_any)
2474 continue;
2476 ineq = get_row_parameter_ineq(tab, row);
2477 if (!ineq)
2478 return -1;
2479 ok = isl_tab_add_ineq(context_tab, ineq->el) >= 0;
2480 isl_vec_free(ineq);
2481 if (!ok)
2482 return -1;
2483 var = &context_tab->con[context_tab->n_con - 1];
2484 if (!context_tab->empty &&
2485 !isl_tab_min_at_most_neg_one(context_tab, var))
2486 r++;
2487 if (isl_tab_rollback(context_tab, snap2) < 0)
2488 return -1;
2490 if (best == -1 || r > best_r) {
2491 best = split;
2492 best_r = r;
2494 if (isl_tab_rollback(context_tab, snap) < 0)
2495 return -1;
2498 return best;
2501 static struct isl_basic_set *context_lex_peek_basic_set(
2502 struct isl_context *context)
2504 struct isl_context_lex *clex = (struct isl_context_lex *)context;
2505 if (!clex->tab)
2506 return NULL;
2507 return isl_tab_peek_bset(clex->tab);
2510 static struct isl_tab *context_lex_peek_tab(struct isl_context *context)
2512 struct isl_context_lex *clex = (struct isl_context_lex *)context;
2513 return clex->tab;
2516 static void context_lex_add_eq(struct isl_context *context, isl_int *eq,
2517 int check, int update)
2519 struct isl_context_lex *clex = (struct isl_context_lex *)context;
2520 if (isl_tab_extend_cons(clex->tab, 2) < 0)
2521 goto error;
2522 if (add_lexmin_eq(clex->tab, eq) < 0)
2523 goto error;
2524 if (check) {
2525 int v = tab_has_valid_sample(clex->tab, eq, 1);
2526 if (v < 0)
2527 goto error;
2528 if (!v)
2529 clex->tab = check_integer_feasible(clex->tab);
2531 if (update)
2532 clex->tab = check_samples(clex->tab, eq, 1);
2533 return;
2534 error:
2535 isl_tab_free(clex->tab);
2536 clex->tab = NULL;
2539 static void context_lex_add_ineq(struct isl_context *context, isl_int *ineq,
2540 int check, int update)
2542 struct isl_context_lex *clex = (struct isl_context_lex *)context;
2543 if (isl_tab_extend_cons(clex->tab, 1) < 0)
2544 goto error;
2545 clex->tab = add_lexmin_ineq(clex->tab, ineq);
2546 if (check) {
2547 int v = tab_has_valid_sample(clex->tab, ineq, 0);
2548 if (v < 0)
2549 goto error;
2550 if (!v)
2551 clex->tab = check_integer_feasible(clex->tab);
2553 if (update)
2554 clex->tab = check_samples(clex->tab, ineq, 0);
2555 return;
2556 error:
2557 isl_tab_free(clex->tab);
2558 clex->tab = NULL;
2561 static isl_stat context_lex_add_ineq_wrap(void *user, isl_int *ineq)
2563 struct isl_context *context = (struct isl_context *)user;
2564 context_lex_add_ineq(context, ineq, 0, 0);
2565 return context->op->is_ok(context) ? isl_stat_ok : isl_stat_error;
2568 /* Check which signs can be obtained by "ineq" on all the currently
2569 * active sample values. See row_sign for more information.
2571 static enum isl_tab_row_sign tab_ineq_sign(struct isl_tab *tab, isl_int *ineq,
2572 int strict)
2574 int i;
2575 int sgn;
2576 isl_int tmp;
2577 enum isl_tab_row_sign res = isl_tab_row_unknown;
2579 isl_assert(tab->mat->ctx, tab->samples, return isl_tab_row_unknown);
2580 isl_assert(tab->mat->ctx, tab->samples->n_col == 1 + tab->n_var,
2581 return isl_tab_row_unknown);
2583 isl_int_init(tmp);
2584 for (i = tab->n_outside; i < tab->n_sample; ++i) {
2585 isl_seq_inner_product(tab->samples->row[i], ineq,
2586 1 + tab->n_var, &tmp);
2587 sgn = isl_int_sgn(tmp);
2588 if (sgn > 0 || (sgn == 0 && strict)) {
2589 if (res == isl_tab_row_unknown)
2590 res = isl_tab_row_pos;
2591 if (res == isl_tab_row_neg)
2592 res = isl_tab_row_any;
2594 if (sgn < 0) {
2595 if (res == isl_tab_row_unknown)
2596 res = isl_tab_row_neg;
2597 if (res == isl_tab_row_pos)
2598 res = isl_tab_row_any;
2600 if (res == isl_tab_row_any)
2601 break;
2603 isl_int_clear(tmp);
2605 return res;
2608 static enum isl_tab_row_sign context_lex_ineq_sign(struct isl_context *context,
2609 isl_int *ineq, int strict)
2611 struct isl_context_lex *clex = (struct isl_context_lex *)context;
2612 return tab_ineq_sign(clex->tab, ineq, strict);
2615 /* Check whether "ineq" can be added to the tableau without rendering
2616 * it infeasible.
2618 static int context_lex_test_ineq(struct isl_context *context, isl_int *ineq)
2620 struct isl_context_lex *clex = (struct isl_context_lex *)context;
2621 struct isl_tab_undo *snap;
2622 int feasible;
2624 if (!clex->tab)
2625 return -1;
2627 if (isl_tab_extend_cons(clex->tab, 1) < 0)
2628 return -1;
2630 snap = isl_tab_snap(clex->tab);
2631 if (isl_tab_push_basis(clex->tab) < 0)
2632 return -1;
2633 clex->tab = add_lexmin_ineq(clex->tab, ineq);
2634 clex->tab = check_integer_feasible(clex->tab);
2635 if (!clex->tab)
2636 return -1;
2637 feasible = !clex->tab->empty;
2638 if (isl_tab_rollback(clex->tab, snap) < 0)
2639 return -1;
2641 return feasible;
2644 static int context_lex_get_div(struct isl_context *context, struct isl_tab *tab,
2645 struct isl_vec *div)
2647 return get_div(tab, context, div);
2650 /* Insert a div specified by "div" to the context tableau at position "pos" and
2651 * return isl_bool_true if the div is obviously non-negative.
2652 * context_tab_add_div will always return isl_bool_true, because all variables
2653 * in a isl_context_lex tableau are non-negative.
2654 * However, if we are using a big parameter in the context, then this only
2655 * reflects the non-negativity of the variable used to _encode_ the
2656 * div, i.e., div' = M + div, so we can't draw any conclusions.
2658 static isl_bool context_lex_insert_div(struct isl_context *context, int pos,
2659 __isl_keep isl_vec *div)
2661 struct isl_context_lex *clex = (struct isl_context_lex *)context;
2662 isl_bool nonneg;
2663 nonneg = context_tab_insert_div(clex->tab, pos, div,
2664 context_lex_add_ineq_wrap, context);
2665 if (nonneg < 0)
2666 return isl_bool_error;
2667 if (clex->tab->M)
2668 return isl_bool_false;
2669 return nonneg;
2672 static int context_lex_detect_equalities(struct isl_context *context,
2673 struct isl_tab *tab)
2675 return 0;
2678 static int context_lex_best_split(struct isl_context *context,
2679 struct isl_tab *tab)
2681 struct isl_context_lex *clex = (struct isl_context_lex *)context;
2682 struct isl_tab_undo *snap;
2683 int r;
2685 snap = isl_tab_snap(clex->tab);
2686 if (isl_tab_push_basis(clex->tab) < 0)
2687 return -1;
2688 r = best_split(tab, clex->tab);
2690 if (r >= 0 && isl_tab_rollback(clex->tab, snap) < 0)
2691 return -1;
2693 return r;
2696 static int context_lex_is_empty(struct isl_context *context)
2698 struct isl_context_lex *clex = (struct isl_context_lex *)context;
2699 if (!clex->tab)
2700 return -1;
2701 return clex->tab->empty;
2704 static void *context_lex_save(struct isl_context *context)
2706 struct isl_context_lex *clex = (struct isl_context_lex *)context;
2707 struct isl_tab_undo *snap;
2709 snap = isl_tab_snap(clex->tab);
2710 if (isl_tab_push_basis(clex->tab) < 0)
2711 return NULL;
2712 if (isl_tab_save_samples(clex->tab) < 0)
2713 return NULL;
2715 return snap;
2718 static void context_lex_restore(struct isl_context *context, void *save)
2720 struct isl_context_lex *clex = (struct isl_context_lex *)context;
2721 if (isl_tab_rollback(clex->tab, (struct isl_tab_undo *)save) < 0) {
2722 isl_tab_free(clex->tab);
2723 clex->tab = NULL;
2727 static void context_lex_discard(void *save)
2731 static int context_lex_is_ok(struct isl_context *context)
2733 struct isl_context_lex *clex = (struct isl_context_lex *)context;
2734 return !!clex->tab;
2737 /* For each variable in the context tableau, check if the variable can
2738 * only attain non-negative values. If so, mark the parameter as non-negative
2739 * in the main tableau. This allows for a more direct identification of some
2740 * cases of violated constraints.
2742 static struct isl_tab *tab_detect_nonnegative_parameters(struct isl_tab *tab,
2743 struct isl_tab *context_tab)
2745 int i;
2746 struct isl_tab_undo *snap;
2747 struct isl_vec *ineq = NULL;
2748 struct isl_tab_var *var;
2749 int n;
2751 if (context_tab->n_var == 0)
2752 return tab;
2754 ineq = isl_vec_alloc(tab->mat->ctx, 1 + context_tab->n_var);
2755 if (!ineq)
2756 goto error;
2758 if (isl_tab_extend_cons(context_tab, 1) < 0)
2759 goto error;
2761 snap = isl_tab_snap(context_tab);
2763 n = 0;
2764 isl_seq_clr(ineq->el, ineq->size);
2765 for (i = 0; i < context_tab->n_var; ++i) {
2766 isl_int_set_si(ineq->el[1 + i], 1);
2767 if (isl_tab_add_ineq(context_tab, ineq->el) < 0)
2768 goto error;
2769 var = &context_tab->con[context_tab->n_con - 1];
2770 if (!context_tab->empty &&
2771 !isl_tab_min_at_most_neg_one(context_tab, var)) {
2772 int j = i;
2773 if (i >= tab->n_param)
2774 j = i - tab->n_param + tab->n_var - tab->n_div;
2775 tab->var[j].is_nonneg = 1;
2776 n++;
2778 isl_int_set_si(ineq->el[1 + i], 0);
2779 if (isl_tab_rollback(context_tab, snap) < 0)
2780 goto error;
2783 if (context_tab->M && n == context_tab->n_var) {
2784 context_tab->mat = isl_mat_drop_cols(context_tab->mat, 2, 1);
2785 context_tab->M = 0;
2788 isl_vec_free(ineq);
2789 return tab;
2790 error:
2791 isl_vec_free(ineq);
2792 isl_tab_free(tab);
2793 return NULL;
2796 static struct isl_tab *context_lex_detect_nonnegative_parameters(
2797 struct isl_context *context, struct isl_tab *tab)
2799 struct isl_context_lex *clex = (struct isl_context_lex *)context;
2800 struct isl_tab_undo *snap;
2802 if (!tab)
2803 return NULL;
2805 snap = isl_tab_snap(clex->tab);
2806 if (isl_tab_push_basis(clex->tab) < 0)
2807 goto error;
2809 tab = tab_detect_nonnegative_parameters(tab, clex->tab);
2811 if (isl_tab_rollback(clex->tab, snap) < 0)
2812 goto error;
2814 return tab;
2815 error:
2816 isl_tab_free(tab);
2817 return NULL;
2820 static void context_lex_invalidate(struct isl_context *context)
2822 struct isl_context_lex *clex = (struct isl_context_lex *)context;
2823 isl_tab_free(clex->tab);
2824 clex->tab = NULL;
2827 static __isl_null struct isl_context *context_lex_free(
2828 struct isl_context *context)
2830 struct isl_context_lex *clex = (struct isl_context_lex *)context;
2831 isl_tab_free(clex->tab);
2832 free(clex);
2834 return NULL;
2837 struct isl_context_op isl_context_lex_op = {
2838 context_lex_detect_nonnegative_parameters,
2839 context_lex_peek_basic_set,
2840 context_lex_peek_tab,
2841 context_lex_add_eq,
2842 context_lex_add_ineq,
2843 context_lex_ineq_sign,
2844 context_lex_test_ineq,
2845 context_lex_get_div,
2846 context_lex_insert_div,
2847 context_lex_detect_equalities,
2848 context_lex_best_split,
2849 context_lex_is_empty,
2850 context_lex_is_ok,
2851 context_lex_save,
2852 context_lex_restore,
2853 context_lex_discard,
2854 context_lex_invalidate,
2855 context_lex_free,
2858 static struct isl_tab *context_tab_for_lexmin(__isl_take isl_basic_set *bset)
2860 struct isl_tab *tab;
2862 if (!bset)
2863 return NULL;
2864 tab = tab_for_lexmin(bset_to_bmap(bset), NULL, 1, 0);
2865 if (isl_tab_track_bset(tab, bset) < 0)
2866 goto error;
2867 tab = isl_tab_init_samples(tab);
2868 return tab;
2869 error:
2870 isl_tab_free(tab);
2871 return NULL;
2874 static struct isl_context *isl_context_lex_alloc(struct isl_basic_set *dom)
2876 struct isl_context_lex *clex;
2878 if (!dom)
2879 return NULL;
2881 clex = isl_alloc_type(dom->ctx, struct isl_context_lex);
2882 if (!clex)
2883 return NULL;
2885 clex->context.op = &isl_context_lex_op;
2887 clex->tab = context_tab_for_lexmin(isl_basic_set_copy(dom));
2888 if (restore_lexmin(clex->tab) < 0)
2889 goto error;
2890 clex->tab = check_integer_feasible(clex->tab);
2891 if (!clex->tab)
2892 goto error;
2894 return &clex->context;
2895 error:
2896 clex->context.op->free(&clex->context);
2897 return NULL;
2900 /* Representation of the context when using generalized basis reduction.
2902 * "shifted" contains the offsets of the unit hypercubes that lie inside the
2903 * context. Any rational point in "shifted" can therefore be rounded
2904 * up to an integer point in the context.
2905 * If the context is constrained by any equality, then "shifted" is not used
2906 * as it would be empty.
2908 struct isl_context_gbr {
2909 struct isl_context context;
2910 struct isl_tab *tab;
2911 struct isl_tab *shifted;
2912 struct isl_tab *cone;
2915 static struct isl_tab *context_gbr_detect_nonnegative_parameters(
2916 struct isl_context *context, struct isl_tab *tab)
2918 struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
2919 if (!tab)
2920 return NULL;
2921 return tab_detect_nonnegative_parameters(tab, cgbr->tab);
2924 static struct isl_basic_set *context_gbr_peek_basic_set(
2925 struct isl_context *context)
2927 struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
2928 if (!cgbr->tab)
2929 return NULL;
2930 return isl_tab_peek_bset(cgbr->tab);
2933 static struct isl_tab *context_gbr_peek_tab(struct isl_context *context)
2935 struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
2936 return cgbr->tab;
2939 /* Initialize the "shifted" tableau of the context, which
2940 * contains the constraints of the original tableau shifted
2941 * by the sum of all negative coefficients. This ensures
2942 * that any rational point in the shifted tableau can
2943 * be rounded up to yield an integer point in the original tableau.
2945 static void gbr_init_shifted(struct isl_context_gbr *cgbr)
2947 int i, j;
2948 struct isl_vec *cst;
2949 struct isl_basic_set *bset = isl_tab_peek_bset(cgbr->tab);
2950 isl_size dim = isl_basic_set_dim(bset, isl_dim_all);
2952 if (dim < 0)
2953 return;
2954 cst = isl_vec_alloc(cgbr->tab->mat->ctx, bset->n_ineq);
2955 if (!cst)
2956 return;
2958 for (i = 0; i < bset->n_ineq; ++i) {
2959 isl_int_set(cst->el[i], bset->ineq[i][0]);
2960 for (j = 0; j < dim; ++j) {
2961 if (!isl_int_is_neg(bset->ineq[i][1 + j]))
2962 continue;
2963 isl_int_add(bset->ineq[i][0], bset->ineq[i][0],
2964 bset->ineq[i][1 + j]);
2968 cgbr->shifted = isl_tab_from_basic_set(bset, 0);
2970 for (i = 0; i < bset->n_ineq; ++i)
2971 isl_int_set(bset->ineq[i][0], cst->el[i]);
2973 isl_vec_free(cst);
2976 /* Check if the shifted tableau is non-empty, and if so
2977 * use the sample point to construct an integer point
2978 * of the context tableau.
2980 static struct isl_vec *gbr_get_shifted_sample(struct isl_context_gbr *cgbr)
2982 struct isl_vec *sample;
2984 if (!cgbr->shifted)
2985 gbr_init_shifted(cgbr);
2986 if (!cgbr->shifted)
2987 return NULL;
2988 if (cgbr->shifted->empty)
2989 return isl_vec_alloc(cgbr->tab->mat->ctx, 0);
2991 sample = isl_tab_get_sample_value(cgbr->shifted);
2992 sample = isl_vec_ceil(sample);
2994 return sample;
2997 static __isl_give isl_basic_set *drop_constant_terms(
2998 __isl_take isl_basic_set *bset)
3000 int i;
3002 if (!bset)
3003 return NULL;
3005 for (i = 0; i < bset->n_eq; ++i)
3006 isl_int_set_si(bset->eq[i][0], 0);
3008 for (i = 0; i < bset->n_ineq; ++i)
3009 isl_int_set_si(bset->ineq[i][0], 0);
3011 return bset;
3014 static int use_shifted(struct isl_context_gbr *cgbr)
3016 if (!cgbr->tab)
3017 return 0;
3018 return cgbr->tab->bmap->n_eq == 0 && cgbr->tab->bmap->n_div == 0;
3021 static struct isl_vec *gbr_get_sample(struct isl_context_gbr *cgbr)
3023 struct isl_basic_set *bset;
3024 struct isl_basic_set *cone;
3026 if (isl_tab_sample_is_integer(cgbr->tab))
3027 return isl_tab_get_sample_value(cgbr->tab);
3029 if (use_shifted(cgbr)) {
3030 struct isl_vec *sample;
3032 sample = gbr_get_shifted_sample(cgbr);
3033 if (!sample || sample->size > 0)
3034 return sample;
3036 isl_vec_free(sample);
3039 if (!cgbr->cone) {
3040 bset = isl_tab_peek_bset(cgbr->tab);
3041 cgbr->cone = isl_tab_from_recession_cone(bset, 0);
3042 if (!cgbr->cone)
3043 return NULL;
3044 if (isl_tab_track_bset(cgbr->cone,
3045 isl_basic_set_copy(bset)) < 0)
3046 return NULL;
3048 if (isl_tab_detect_implicit_equalities(cgbr->cone) < 0)
3049 return NULL;
3051 if (cgbr->cone->n_dead == cgbr->cone->n_col) {
3052 struct isl_vec *sample;
3053 struct isl_tab_undo *snap;
3055 if (cgbr->tab->basis) {
3056 if (cgbr->tab->basis->n_col != 1 + cgbr->tab->n_var) {
3057 isl_mat_free(cgbr->tab->basis);
3058 cgbr->tab->basis = NULL;
3060 cgbr->tab->n_zero = 0;
3061 cgbr->tab->n_unbounded = 0;
3064 snap = isl_tab_snap(cgbr->tab);
3066 sample = isl_tab_sample(cgbr->tab);
3068 if (!sample || isl_tab_rollback(cgbr->tab, snap) < 0) {
3069 isl_vec_free(sample);
3070 return NULL;
3073 return sample;
3076 cone = isl_basic_set_dup(isl_tab_peek_bset(cgbr->cone));
3077 cone = drop_constant_terms(cone);
3078 cone = isl_basic_set_update_from_tab(cone, cgbr->cone);
3079 cone = isl_basic_set_underlying_set(cone);
3080 cone = isl_basic_set_gauss(cone, NULL);
3082 bset = isl_basic_set_dup(isl_tab_peek_bset(cgbr->tab));
3083 bset = isl_basic_set_update_from_tab(bset, cgbr->tab);
3084 bset = isl_basic_set_underlying_set(bset);
3085 bset = isl_basic_set_gauss(bset, NULL);
3087 return isl_basic_set_sample_with_cone(bset, cone);
3090 static void check_gbr_integer_feasible(struct isl_context_gbr *cgbr)
3092 struct isl_vec *sample;
3094 if (!cgbr->tab)
3095 return;
3097 if (cgbr->tab->empty)
3098 return;
3100 sample = gbr_get_sample(cgbr);
3101 if (!sample)
3102 goto error;
3104 if (sample->size == 0) {
3105 isl_vec_free(sample);
3106 if (isl_tab_mark_empty(cgbr->tab) < 0)
3107 goto error;
3108 return;
3111 if (isl_tab_add_sample(cgbr->tab, sample) < 0)
3112 goto error;
3114 return;
3115 error:
3116 isl_tab_free(cgbr->tab);
3117 cgbr->tab = NULL;
3120 static struct isl_tab *add_gbr_eq(struct isl_tab *tab, isl_int *eq)
3122 if (!tab)
3123 return NULL;
3125 if (isl_tab_extend_cons(tab, 2) < 0)
3126 goto error;
3128 if (isl_tab_add_eq(tab, eq) < 0)
3129 goto error;
3131 return tab;
3132 error:
3133 isl_tab_free(tab);
3134 return NULL;
3137 /* Add the equality described by "eq" to the context.
3138 * If "check" is set, then we check if the context is empty after
3139 * adding the equality.
3140 * If "update" is set, then we check if the samples are still valid.
3142 * We do not explicitly add shifted copies of the equality to
3143 * cgbr->shifted since they would conflict with each other.
3144 * Instead, we directly mark cgbr->shifted empty.
3146 static void context_gbr_add_eq(struct isl_context *context, isl_int *eq,
3147 int check, int update)
3149 struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3151 cgbr->tab = add_gbr_eq(cgbr->tab, eq);
3153 if (cgbr->shifted && !cgbr->shifted->empty && use_shifted(cgbr)) {
3154 if (isl_tab_mark_empty(cgbr->shifted) < 0)
3155 goto error;
3158 if (cgbr->cone && cgbr->cone->n_col != cgbr->cone->n_dead) {
3159 if (isl_tab_extend_cons(cgbr->cone, 2) < 0)
3160 goto error;
3161 if (isl_tab_add_eq(cgbr->cone, eq) < 0)
3162 goto error;
3165 if (check) {
3166 int v = tab_has_valid_sample(cgbr->tab, eq, 1);
3167 if (v < 0)
3168 goto error;
3169 if (!v)
3170 check_gbr_integer_feasible(cgbr);
3172 if (update)
3173 cgbr->tab = check_samples(cgbr->tab, eq, 1);
3174 return;
3175 error:
3176 isl_tab_free(cgbr->tab);
3177 cgbr->tab = NULL;
3180 static void add_gbr_ineq(struct isl_context_gbr *cgbr, isl_int *ineq)
3182 if (!cgbr->tab)
3183 return;
3185 if (isl_tab_extend_cons(cgbr->tab, 1) < 0)
3186 goto error;
3188 if (isl_tab_add_ineq(cgbr->tab, ineq) < 0)
3189 goto error;
3191 if (cgbr->shifted && !cgbr->shifted->empty && use_shifted(cgbr)) {
3192 int i;
3193 isl_size dim;
3194 dim = isl_basic_map_dim(cgbr->tab->bmap, isl_dim_all);
3195 if (dim < 0)
3196 goto error;
3198 if (isl_tab_extend_cons(cgbr->shifted, 1) < 0)
3199 goto error;
3201 for (i = 0; i < dim; ++i) {
3202 if (!isl_int_is_neg(ineq[1 + i]))
3203 continue;
3204 isl_int_add(ineq[0], ineq[0], ineq[1 + i]);
3207 if (isl_tab_add_ineq(cgbr->shifted, ineq) < 0)
3208 goto error;
3210 for (i = 0; i < dim; ++i) {
3211 if (!isl_int_is_neg(ineq[1 + i]))
3212 continue;
3213 isl_int_sub(ineq[0], ineq[0], ineq[1 + i]);
3217 if (cgbr->cone && cgbr->cone->n_col != cgbr->cone->n_dead) {
3218 if (isl_tab_extend_cons(cgbr->cone, 1) < 0)
3219 goto error;
3220 if (isl_tab_add_ineq(cgbr->cone, ineq) < 0)
3221 goto error;
3224 return;
3225 error:
3226 isl_tab_free(cgbr->tab);
3227 cgbr->tab = NULL;
3230 static void context_gbr_add_ineq(struct isl_context *context, isl_int *ineq,
3231 int check, int update)
3233 struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3235 add_gbr_ineq(cgbr, ineq);
3236 if (!cgbr->tab)
3237 return;
3239 if (check) {
3240 int v = tab_has_valid_sample(cgbr->tab, ineq, 0);
3241 if (v < 0)
3242 goto error;
3243 if (!v)
3244 check_gbr_integer_feasible(cgbr);
3246 if (update)
3247 cgbr->tab = check_samples(cgbr->tab, ineq, 0);
3248 return;
3249 error:
3250 isl_tab_free(cgbr->tab);
3251 cgbr->tab = NULL;
3254 static isl_stat context_gbr_add_ineq_wrap(void *user, isl_int *ineq)
3256 struct isl_context *context = (struct isl_context *)user;
3257 context_gbr_add_ineq(context, ineq, 0, 0);
3258 return context->op->is_ok(context) ? isl_stat_ok : isl_stat_error;
3261 static enum isl_tab_row_sign context_gbr_ineq_sign(struct isl_context *context,
3262 isl_int *ineq, int strict)
3264 struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3265 return tab_ineq_sign(cgbr->tab, ineq, strict);
3268 /* Check whether "ineq" can be added to the tableau without rendering
3269 * it infeasible.
3271 static int context_gbr_test_ineq(struct isl_context *context, isl_int *ineq)
3273 struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3274 struct isl_tab_undo *snap;
3275 struct isl_tab_undo *shifted_snap = NULL;
3276 struct isl_tab_undo *cone_snap = NULL;
3277 int feasible;
3279 if (!cgbr->tab)
3280 return -1;
3282 if (isl_tab_extend_cons(cgbr->tab, 1) < 0)
3283 return -1;
3285 snap = isl_tab_snap(cgbr->tab);
3286 if (cgbr->shifted)
3287 shifted_snap = isl_tab_snap(cgbr->shifted);
3288 if (cgbr->cone)
3289 cone_snap = isl_tab_snap(cgbr->cone);
3290 add_gbr_ineq(cgbr, ineq);
3291 check_gbr_integer_feasible(cgbr);
3292 if (!cgbr->tab)
3293 return -1;
3294 feasible = !cgbr->tab->empty;
3295 if (isl_tab_rollback(cgbr->tab, snap) < 0)
3296 return -1;
3297 if (shifted_snap) {
3298 if (isl_tab_rollback(cgbr->shifted, shifted_snap))
3299 return -1;
3300 } else if (cgbr->shifted) {
3301 isl_tab_free(cgbr->shifted);
3302 cgbr->shifted = NULL;
3304 if (cone_snap) {
3305 if (isl_tab_rollback(cgbr->cone, cone_snap))
3306 return -1;
3307 } else if (cgbr->cone) {
3308 isl_tab_free(cgbr->cone);
3309 cgbr->cone = NULL;
3312 return feasible;
3315 /* Return the column of the last of the variables associated to
3316 * a column that has a non-zero coefficient.
3317 * This function is called in a context where only coefficients
3318 * of parameters or divs can be non-zero.
3320 static int last_non_zero_var_col(struct isl_tab *tab, isl_int *p)
3322 int i;
3323 int col;
3325 if (tab->n_var == 0)
3326 return -1;
3328 for (i = tab->n_var - 1; i >= 0; --i) {
3329 if (i >= tab->n_param && i < tab->n_var - tab->n_div)
3330 continue;
3331 if (tab->var[i].is_row)
3332 continue;
3333 col = tab->var[i].index;
3334 if (!isl_int_is_zero(p[col]))
3335 return col;
3338 return -1;
3341 /* Look through all the recently added equalities in the context
3342 * to see if we can propagate any of them to the main tableau.
3344 * The newly added equalities in the context are encoded as pairs
3345 * of inequalities starting at inequality "first".
3347 * We tentatively add each of these equalities to the main tableau
3348 * and if this happens to result in a row with a final coefficient
3349 * that is one or negative one, we use it to kill a column
3350 * in the main tableau. Otherwise, we discard the tentatively
3351 * added row.
3352 * This tentative addition of equality constraints turns
3353 * on the undo facility of the tableau. Turn it off again
3354 * at the end, assuming it was turned off to begin with.
3356 * Return 0 on success and -1 on failure.
3358 static int propagate_equalities(struct isl_context_gbr *cgbr,
3359 struct isl_tab *tab, unsigned first)
3361 int i;
3362 struct isl_vec *eq = NULL;
3363 isl_bool needs_undo;
3365 needs_undo = isl_tab_need_undo(tab);
3366 if (needs_undo < 0)
3367 goto error;
3368 eq = isl_vec_alloc(tab->mat->ctx, 1 + tab->n_var);
3369 if (!eq)
3370 goto error;
3372 if (isl_tab_extend_cons(tab, (cgbr->tab->bmap->n_ineq - first)/2) < 0)
3373 goto error;
3375 isl_seq_clr(eq->el + 1 + tab->n_param,
3376 tab->n_var - tab->n_param - tab->n_div);
3377 for (i = first; i < cgbr->tab->bmap->n_ineq; i += 2) {
3378 int j;
3379 int r;
3380 struct isl_tab_undo *snap;
3381 snap = isl_tab_snap(tab);
3383 isl_seq_cpy(eq->el, cgbr->tab->bmap->ineq[i], 1 + tab->n_param);
3384 isl_seq_cpy(eq->el + 1 + tab->n_var - tab->n_div,
3385 cgbr->tab->bmap->ineq[i] + 1 + tab->n_param,
3386 tab->n_div);
3388 r = isl_tab_add_row(tab, eq->el);
3389 if (r < 0)
3390 goto error;
3391 r = tab->con[r].index;
3392 j = last_non_zero_var_col(tab, tab->mat->row[r] + 2 + tab->M);
3393 if (j < 0 || j < tab->n_dead ||
3394 !isl_int_is_one(tab->mat->row[r][0]) ||
3395 (!isl_int_is_one(tab->mat->row[r][2 + tab->M + j]) &&
3396 !isl_int_is_negone(tab->mat->row[r][2 + tab->M + j]))) {
3397 if (isl_tab_rollback(tab, snap) < 0)
3398 goto error;
3399 continue;
3401 if (isl_tab_pivot(tab, r, j) < 0)
3402 goto error;
3403 if (isl_tab_kill_col(tab, j) < 0)
3404 goto error;
3406 if (restore_lexmin(tab) < 0)
3407 goto error;
3410 if (!needs_undo)
3411 isl_tab_clear_undo(tab);
3412 isl_vec_free(eq);
3414 return 0;
3415 error:
3416 isl_vec_free(eq);
3417 isl_tab_free(cgbr->tab);
3418 cgbr->tab = NULL;
3419 return -1;
3422 static int context_gbr_detect_equalities(struct isl_context *context,
3423 struct isl_tab *tab)
3425 struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3426 unsigned n_ineq;
3428 if (!cgbr->cone) {
3429 struct isl_basic_set *bset = isl_tab_peek_bset(cgbr->tab);
3430 cgbr->cone = isl_tab_from_recession_cone(bset, 0);
3431 if (!cgbr->cone)
3432 goto error;
3433 if (isl_tab_track_bset(cgbr->cone,
3434 isl_basic_set_copy(bset)) < 0)
3435 goto error;
3437 if (isl_tab_detect_implicit_equalities(cgbr->cone) < 0)
3438 goto error;
3440 n_ineq = cgbr->tab->bmap->n_ineq;
3441 cgbr->tab = isl_tab_detect_equalities(cgbr->tab, cgbr->cone);
3442 if (!cgbr->tab)
3443 return -1;
3444 if (cgbr->tab->bmap->n_ineq > n_ineq &&
3445 propagate_equalities(cgbr, tab, n_ineq) < 0)
3446 return -1;
3448 return 0;
3449 error:
3450 isl_tab_free(cgbr->tab);
3451 cgbr->tab = NULL;
3452 return -1;
3455 static int context_gbr_get_div(struct isl_context *context, struct isl_tab *tab,
3456 struct isl_vec *div)
3458 return get_div(tab, context, div);
3461 static isl_bool context_gbr_insert_div(struct isl_context *context, int pos,
3462 __isl_keep isl_vec *div)
3464 struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3465 if (cgbr->cone) {
3466 int r, o_div;
3467 isl_size n_div;
3469 n_div = isl_basic_map_dim(cgbr->cone->bmap, isl_dim_div);
3470 if (n_div < 0)
3471 return isl_bool_error;
3472 o_div = cgbr->cone->n_var - n_div;
3474 if (isl_tab_extend_cons(cgbr->cone, 3) < 0)
3475 return isl_bool_error;
3476 if (isl_tab_extend_vars(cgbr->cone, 1) < 0)
3477 return isl_bool_error;
3478 if ((r = isl_tab_insert_var(cgbr->cone, pos)) <0)
3479 return isl_bool_error;
3481 cgbr->cone->bmap = isl_basic_map_insert_div(cgbr->cone->bmap,
3482 r - o_div, div);
3483 if (!cgbr->cone->bmap)
3484 return isl_bool_error;
3485 if (isl_tab_push_var(cgbr->cone, isl_tab_undo_bmap_div,
3486 &cgbr->cone->var[r]) < 0)
3487 return isl_bool_error;
3489 return context_tab_insert_div(cgbr->tab, pos, div,
3490 context_gbr_add_ineq_wrap, context);
3493 static int context_gbr_best_split(struct isl_context *context,
3494 struct isl_tab *tab)
3496 struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3497 struct isl_tab_undo *snap;
3498 int r;
3500 snap = isl_tab_snap(cgbr->tab);
3501 r = best_split(tab, cgbr->tab);
3503 if (r >= 0 && isl_tab_rollback(cgbr->tab, snap) < 0)
3504 return -1;
3506 return r;
3509 static int context_gbr_is_empty(struct isl_context *context)
3511 struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3512 if (!cgbr->tab)
3513 return -1;
3514 return cgbr->tab->empty;
3517 struct isl_gbr_tab_undo {
3518 struct isl_tab_undo *tab_snap;
3519 struct isl_tab_undo *shifted_snap;
3520 struct isl_tab_undo *cone_snap;
3523 static void *context_gbr_save(struct isl_context *context)
3525 struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3526 struct isl_gbr_tab_undo *snap;
3528 if (!cgbr->tab)
3529 return NULL;
3531 snap = isl_alloc_type(cgbr->tab->mat->ctx, struct isl_gbr_tab_undo);
3532 if (!snap)
3533 return NULL;
3535 snap->tab_snap = isl_tab_snap(cgbr->tab);
3536 if (isl_tab_save_samples(cgbr->tab) < 0)
3537 goto error;
3539 if (cgbr->shifted)
3540 snap->shifted_snap = isl_tab_snap(cgbr->shifted);
3541 else
3542 snap->shifted_snap = NULL;
3544 if (cgbr->cone)
3545 snap->cone_snap = isl_tab_snap(cgbr->cone);
3546 else
3547 snap->cone_snap = NULL;
3549 return snap;
3550 error:
3551 free(snap);
3552 return NULL;
3555 static void context_gbr_restore(struct isl_context *context, void *save)
3557 struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3558 struct isl_gbr_tab_undo *snap = (struct isl_gbr_tab_undo *)save;
3559 if (!snap)
3560 goto error;
3561 if (isl_tab_rollback(cgbr->tab, snap->tab_snap) < 0)
3562 goto error;
3564 if (snap->shifted_snap) {
3565 if (isl_tab_rollback(cgbr->shifted, snap->shifted_snap) < 0)
3566 goto error;
3567 } else if (cgbr->shifted) {
3568 isl_tab_free(cgbr->shifted);
3569 cgbr->shifted = NULL;
3572 if (snap->cone_snap) {
3573 if (isl_tab_rollback(cgbr->cone, snap->cone_snap) < 0)
3574 goto error;
3575 } else if (cgbr->cone) {
3576 isl_tab_free(cgbr->cone);
3577 cgbr->cone = NULL;
3580 free(snap);
3582 return;
3583 error:
3584 free(snap);
3585 isl_tab_free(cgbr->tab);
3586 cgbr->tab = NULL;
3589 static void context_gbr_discard(void *save)
3591 struct isl_gbr_tab_undo *snap = (struct isl_gbr_tab_undo *)save;
3592 free(snap);
3595 static int context_gbr_is_ok(struct isl_context *context)
3597 struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3598 return !!cgbr->tab;
3601 static void context_gbr_invalidate(struct isl_context *context)
3603 struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3604 isl_tab_free(cgbr->tab);
3605 cgbr->tab = NULL;
3608 static __isl_null struct isl_context *context_gbr_free(
3609 struct isl_context *context)
3611 struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3612 isl_tab_free(cgbr->tab);
3613 isl_tab_free(cgbr->shifted);
3614 isl_tab_free(cgbr->cone);
3615 free(cgbr);
3617 return NULL;
3620 struct isl_context_op isl_context_gbr_op = {
3621 context_gbr_detect_nonnegative_parameters,
3622 context_gbr_peek_basic_set,
3623 context_gbr_peek_tab,
3624 context_gbr_add_eq,
3625 context_gbr_add_ineq,
3626 context_gbr_ineq_sign,
3627 context_gbr_test_ineq,
3628 context_gbr_get_div,
3629 context_gbr_insert_div,
3630 context_gbr_detect_equalities,
3631 context_gbr_best_split,
3632 context_gbr_is_empty,
3633 context_gbr_is_ok,
3634 context_gbr_save,
3635 context_gbr_restore,
3636 context_gbr_discard,
3637 context_gbr_invalidate,
3638 context_gbr_free,
3641 static struct isl_context *isl_context_gbr_alloc(__isl_keep isl_basic_set *dom)
3643 struct isl_context_gbr *cgbr;
3645 if (!dom)
3646 return NULL;
3648 cgbr = isl_calloc_type(dom->ctx, struct isl_context_gbr);
3649 if (!cgbr)
3650 return NULL;
3652 cgbr->context.op = &isl_context_gbr_op;
3654 cgbr->shifted = NULL;
3655 cgbr->cone = NULL;
3656 cgbr->tab = isl_tab_from_basic_set(dom, 1);
3657 cgbr->tab = isl_tab_init_samples(cgbr->tab);
3658 if (!cgbr->tab)
3659 goto error;
3660 check_gbr_integer_feasible(cgbr);
3662 return &cgbr->context;
3663 error:
3664 cgbr->context.op->free(&cgbr->context);
3665 return NULL;
3668 /* Allocate a context corresponding to "dom".
3669 * The representation specific fields are initialized by
3670 * isl_context_lex_alloc or isl_context_gbr_alloc.
3671 * The shared "n_unknown" field is initialized to the number
3672 * of final unknown integer divisions in "dom".
3674 static struct isl_context *isl_context_alloc(__isl_keep isl_basic_set *dom)
3676 struct isl_context *context;
3677 int first;
3678 isl_size n_div;
3680 if (!dom)
3681 return NULL;
3683 if (dom->ctx->opt->context == ISL_CONTEXT_LEXMIN)
3684 context = isl_context_lex_alloc(dom);
3685 else
3686 context = isl_context_gbr_alloc(dom);
3688 if (!context)
3689 return NULL;
3691 first = isl_basic_set_first_unknown_div(dom);
3692 n_div = isl_basic_set_dim(dom, isl_dim_div);
3693 if (first < 0 || n_div < 0)
3694 return context->op->free(context);
3695 context->n_unknown = n_div - first;
3697 return context;
3700 /* Initialize some common fields of "sol", which keeps track
3701 * of the solution of an optimization problem on "bmap" over
3702 * the domain "dom".
3703 * If "max" is set, then a maximization problem is being solved, rather than
3704 * a minimization problem, which means that the variables in the
3705 * tableau have value "M - x" rather than "M + x".
3707 static isl_stat sol_init(struct isl_sol *sol, __isl_keep isl_basic_map *bmap,
3708 __isl_keep isl_basic_set *dom, int max)
3710 sol->rational = ISL_F_ISSET(bmap, ISL_BASIC_MAP_RATIONAL);
3711 sol->dec_level.callback.run = &sol_dec_level_wrap;
3712 sol->dec_level.sol = sol;
3713 sol->max = max;
3714 sol->n_out = isl_basic_map_dim(bmap, isl_dim_out);
3715 sol->space = isl_basic_map_get_space(bmap);
3717 sol->context = isl_context_alloc(dom);
3718 if (sol->n_out < 0 || !sol->space || !sol->context)
3719 return isl_stat_error;
3721 return isl_stat_ok;
3724 /* Construct an isl_sol_map structure for accumulating the solution.
3725 * If track_empty is set, then we also keep track of the parts
3726 * of the context where there is no solution.
3727 * If max is set, then we are solving a maximization, rather than
3728 * a minimization problem, which means that the variables in the
3729 * tableau have value "M - x" rather than "M + x".
3731 static struct isl_sol *sol_map_init(__isl_keep isl_basic_map *bmap,
3732 __isl_take isl_basic_set *dom, int track_empty, int max)
3734 struct isl_sol_map *sol_map = NULL;
3735 isl_space *space;
3737 if (!bmap)
3738 goto error;
3740 sol_map = isl_calloc_type(bmap->ctx, struct isl_sol_map);
3741 if (!sol_map)
3742 goto error;
3744 sol_map->sol.free = &sol_map_free;
3745 if (sol_init(&sol_map->sol, bmap, dom, max) < 0)
3746 goto error;
3747 sol_map->sol.add = &sol_map_add_wrap;
3748 sol_map->sol.add_empty = track_empty ? &sol_map_add_empty_wrap : NULL;
3749 space = isl_space_copy(sol_map->sol.space);
3750 sol_map->map = isl_map_alloc_space(space, 1, ISL_MAP_DISJOINT);
3751 if (!sol_map->map)
3752 goto error;
3754 if (track_empty) {
3755 sol_map->empty = isl_set_alloc_space(isl_basic_set_get_space(dom),
3756 1, ISL_SET_DISJOINT);
3757 if (!sol_map->empty)
3758 goto error;
3761 isl_basic_set_free(dom);
3762 return &sol_map->sol;
3763 error:
3764 isl_basic_set_free(dom);
3765 sol_free(&sol_map->sol);
3766 return NULL;
3769 /* Check whether all coefficients of (non-parameter) variables
3770 * are non-positive, meaning that no pivots can be performed on the row.
3772 static int is_critical(struct isl_tab *tab, int row)
3774 int j;
3775 unsigned off = 2 + tab->M;
3777 for (j = tab->n_dead; j < tab->n_col; ++j) {
3778 if (col_is_parameter_var(tab, j))
3779 continue;
3781 if (isl_int_is_pos(tab->mat->row[row][off + j]))
3782 return 0;
3785 return 1;
3788 /* Check whether the inequality represented by vec is strict over the integers,
3789 * i.e., there are no integer values satisfying the constraint with
3790 * equality. This happens if the gcd of the coefficients is not a divisor
3791 * of the constant term. If so, scale the constraint down by the gcd
3792 * of the coefficients.
3794 static int is_strict(struct isl_vec *vec)
3796 isl_int gcd;
3797 int strict = 0;
3799 isl_int_init(gcd);
3800 isl_seq_gcd(vec->el + 1, vec->size - 1, &gcd);
3801 if (!isl_int_is_one(gcd)) {
3802 strict = !isl_int_is_divisible_by(vec->el[0], gcd);
3803 isl_int_fdiv_q(vec->el[0], vec->el[0], gcd);
3804 isl_seq_scale_down(vec->el + 1, vec->el + 1, gcd, vec->size-1);
3806 isl_int_clear(gcd);
3808 return strict;
3811 /* Determine the sign of the given row of the main tableau.
3812 * The result is one of
3813 * isl_tab_row_pos: always non-negative; no pivot needed
3814 * isl_tab_row_neg: always non-positive; pivot
3815 * isl_tab_row_any: can be both positive and negative; split
3817 * We first handle some simple cases
3818 * - the row sign may be known already
3819 * - the row may be obviously non-negative
3820 * - the parametric constant may be equal to that of another row
3821 * for which we know the sign. This sign will be either "pos" or
3822 * "any". If it had been "neg" then we would have pivoted before.
3824 * If none of these cases hold, we check the value of the row for each
3825 * of the currently active samples. Based on the signs of these values
3826 * we make an initial determination of the sign of the row.
3828 * all zero -> unk(nown)
3829 * all non-negative -> pos
3830 * all non-positive -> neg
3831 * both negative and positive -> all
3833 * If we end up with "all", we are done.
3834 * Otherwise, we perform a check for positive and/or negative
3835 * values as follows.
3837 * samples neg unk pos
3838 * <0 ? Y N Y N
3839 * pos any pos
3840 * >0 ? Y N Y N
3841 * any neg any neg
3843 * There is no special sign for "zero", because we can usually treat zero
3844 * as either non-negative or non-positive, whatever works out best.
3845 * However, if the row is "critical", meaning that pivoting is impossible
3846 * then we don't want to limp zero with the non-positive case, because
3847 * then we we would lose the solution for those values of the parameters
3848 * where the value of the row is zero. Instead, we treat 0 as non-negative
3849 * ensuring a split if the row can attain both zero and negative values.
3850 * The same happens when the original constraint was one that could not
3851 * be satisfied with equality by any integer values of the parameters.
3852 * In this case, we normalize the constraint, but then a value of zero
3853 * for the normalized constraint is actually a positive value for the
3854 * original constraint, so again we need to treat zero as non-negative.
3855 * In both these cases, we have the following decision tree instead:
3857 * all non-negative -> pos
3858 * all negative -> neg
3859 * both negative and non-negative -> all
3861 * samples neg pos
3862 * <0 ? Y N
3863 * any pos
3864 * >=0 ? Y N
3865 * any neg
3867 static enum isl_tab_row_sign row_sign(struct isl_tab *tab,
3868 struct isl_sol *sol, int row)
3870 struct isl_vec *ineq = NULL;
3871 enum isl_tab_row_sign res = isl_tab_row_unknown;
3872 int critical;
3873 int strict;
3874 int row2;
3876 if (tab->row_sign[row] != isl_tab_row_unknown)
3877 return tab->row_sign[row];
3878 if (is_obviously_nonneg(tab, row))
3879 return isl_tab_row_pos;
3880 for (row2 = tab->n_redundant; row2 < tab->n_row; ++row2) {
3881 if (tab->row_sign[row2] == isl_tab_row_unknown)
3882 continue;
3883 if (identical_parameter_line(tab, row, row2))
3884 return tab->row_sign[row2];
3887 critical = is_critical(tab, row);
3889 ineq = get_row_parameter_ineq(tab, row);
3890 if (!ineq)
3891 goto error;
3893 strict = is_strict(ineq);
3895 res = sol->context->op->ineq_sign(sol->context, ineq->el,
3896 critical || strict);
3898 if (res == isl_tab_row_unknown || res == isl_tab_row_pos) {
3899 /* test for negative values */
3900 int feasible;
3901 isl_seq_neg(ineq->el, ineq->el, ineq->size);
3902 isl_int_sub_ui(ineq->el[0], ineq->el[0], 1);
3904 feasible = sol->context->op->test_ineq(sol->context, ineq->el);
3905 if (feasible < 0)
3906 goto error;
3907 if (!feasible)
3908 res = isl_tab_row_pos;
3909 else
3910 res = (res == isl_tab_row_unknown) ? isl_tab_row_neg
3911 : isl_tab_row_any;
3912 if (res == isl_tab_row_neg) {
3913 isl_seq_neg(ineq->el, ineq->el, ineq->size);
3914 isl_int_sub_ui(ineq->el[0], ineq->el[0], 1);
3918 if (res == isl_tab_row_neg) {
3919 /* test for positive values */
3920 int feasible;
3921 if (!critical && !strict)
3922 isl_int_sub_ui(ineq->el[0], ineq->el[0], 1);
3924 feasible = sol->context->op->test_ineq(sol->context, ineq->el);
3925 if (feasible < 0)
3926 goto error;
3927 if (feasible)
3928 res = isl_tab_row_any;
3931 isl_vec_free(ineq);
3932 return res;
3933 error:
3934 isl_vec_free(ineq);
3935 return isl_tab_row_unknown;
3938 static void find_solutions(struct isl_sol *sol, struct isl_tab *tab);
3940 /* Find solutions for values of the parameters that satisfy the given
3941 * inequality.
3943 * We currently take a snapshot of the context tableau that is reset
3944 * when we return from this function, while we make a copy of the main
3945 * tableau, leaving the original main tableau untouched.
3946 * These are fairly arbitrary choices. Making a copy also of the context
3947 * tableau would obviate the need to undo any changes made to it later,
3948 * while taking a snapshot of the main tableau could reduce memory usage.
3949 * If we were to switch to taking a snapshot of the main tableau,
3950 * we would have to keep in mind that we need to save the row signs
3951 * and that we need to do this before saving the current basis
3952 * such that the basis has been restore before we restore the row signs.
3954 static void find_in_pos(struct isl_sol *sol, struct isl_tab *tab, isl_int *ineq)
3956 void *saved;
3958 if (!sol->context)
3959 goto error;
3961 tab = isl_tab_dup(tab);
3962 if (!tab)
3963 goto error;
3965 saved = sol->context->op->save(sol->context);
3967 sol_context_add_ineq(sol, ineq, 0, 1);
3969 find_solutions(sol, tab);
3971 if (!sol->error)
3972 sol->context->op->restore(sol->context, saved);
3973 else
3974 sol->context->op->discard(saved);
3975 return;
3976 error:
3977 sol->error = 1;
3980 /* Record the absence of solutions for those values of the parameters
3981 * that do not satisfy the given inequality with equality.
3983 static void no_sol_in_strict(struct isl_sol *sol,
3984 struct isl_tab *tab, struct isl_vec *ineq)
3986 int empty;
3987 void *saved;
3989 if (!sol->context || sol->error)
3990 goto error;
3991 saved = sol->context->op->save(sol->context);
3993 isl_int_sub_ui(ineq->el[0], ineq->el[0], 1);
3995 sol_context_add_ineq(sol, ineq->el, 1, 0);
3997 empty = tab->empty;
3998 tab->empty = 1;
3999 sol_add(sol, tab);
4000 tab->empty = empty;
4002 isl_int_add_ui(ineq->el[0], ineq->el[0], 1);
4004 sol->context->op->restore(sol->context, saved);
4005 if (!sol->context->op->is_ok(sol->context))
4006 goto error;
4007 return;
4008 error:
4009 sol->error = 1;
4012 /* Reset all row variables that are marked to have a sign that may
4013 * be both positive and negative to have an unknown sign.
4015 static void reset_any_to_unknown(struct isl_tab *tab)
4017 int row;
4019 for (row = tab->n_redundant; row < tab->n_row; ++row) {
4020 if (!isl_tab_var_from_row(tab, row)->is_nonneg)
4021 continue;
4022 if (tab->row_sign[row] == isl_tab_row_any)
4023 tab->row_sign[row] = isl_tab_row_unknown;
4027 /* Compute the lexicographic minimum of the set represented by the main
4028 * tableau "tab" within the context "sol->context_tab".
4029 * On entry the sample value of the main tableau is lexicographically
4030 * less than or equal to this lexicographic minimum.
4031 * Pivots are performed until a feasible point is found, which is then
4032 * necessarily equal to the minimum, or until the tableau is found to
4033 * be infeasible. Some pivots may need to be performed for only some
4034 * feasible values of the context tableau. If so, the context tableau
4035 * is split into a part where the pivot is needed and a part where it is not.
4037 * Whenever we enter the main loop, the main tableau is such that no
4038 * "obvious" pivots need to be performed on it, where "obvious" means
4039 * that the given row can be seen to be negative without looking at
4040 * the context tableau. In particular, for non-parametric problems,
4041 * no pivots need to be performed on the main tableau.
4042 * The caller of find_solutions is responsible for making this property
4043 * hold prior to the first iteration of the loop, while restore_lexmin
4044 * is called before every other iteration.
4046 * Inside the main loop, we first examine the signs of the rows of
4047 * the main tableau within the context of the context tableau.
4048 * If we find a row that is always non-positive for all values of
4049 * the parameters satisfying the context tableau and negative for at
4050 * least one value of the parameters, we perform the appropriate pivot
4051 * and start over. An exception is the case where no pivot can be
4052 * performed on the row. In this case, we require that the sign of
4053 * the row is negative for all values of the parameters (rather than just
4054 * non-positive). This special case is handled inside row_sign, which
4055 * will say that the row can have any sign if it determines that it can
4056 * attain both negative and zero values.
4058 * If we can't find a row that always requires a pivot, but we can find
4059 * one or more rows that require a pivot for some values of the parameters
4060 * (i.e., the row can attain both positive and negative signs), then we split
4061 * the context tableau into two parts, one where we force the sign to be
4062 * non-negative and one where we force is to be negative.
4063 * The non-negative part is handled by a recursive call (through find_in_pos).
4064 * Upon returning from this call, we continue with the negative part and
4065 * perform the required pivot.
4067 * If no such rows can be found, all rows are non-negative and we have
4068 * found a (rational) feasible point. If we only wanted a rational point
4069 * then we are done.
4070 * Otherwise, we check if all values of the sample point of the tableau
4071 * are integral for the variables. If so, we have found the minimal
4072 * integral point and we are done.
4073 * If the sample point is not integral, then we need to make a distinction
4074 * based on whether the constant term is non-integral or the coefficients
4075 * of the parameters. Furthermore, in order to decide how to handle
4076 * the non-integrality, we also need to know whether the coefficients
4077 * of the other columns in the tableau are integral. This leads
4078 * to the following table. The first two rows do not correspond
4079 * to a non-integral sample point and are only mentioned for completeness.
4081 * constant parameters other
4083 * int int int |
4084 * int int rat | -> no problem
4086 * rat int int -> fail
4088 * rat int rat -> cut
4090 * int rat rat |
4091 * rat rat rat | -> parametric cut
4093 * int rat int |
4094 * rat rat int | -> split context
4096 * If the parametric constant is completely integral, then there is nothing
4097 * to be done. If the constant term is non-integral, but all the other
4098 * coefficient are integral, then there is nothing that can be done
4099 * and the tableau has no integral solution.
4100 * If, on the other hand, one or more of the other columns have rational
4101 * coefficients, but the parameter coefficients are all integral, then
4102 * we can perform a regular (non-parametric) cut.
4103 * Finally, if there is any parameter coefficient that is non-integral,
4104 * then we need to involve the context tableau. There are two cases here.
4105 * If at least one other column has a rational coefficient, then we
4106 * can perform a parametric cut in the main tableau by adding a new
4107 * integer division in the context tableau.
4108 * If all other columns have integral coefficients, then we need to
4109 * enforce that the rational combination of parameters (c + \sum a_i y_i)/m
4110 * is always integral. We do this by introducing an integer division
4111 * q = floor((c + \sum a_i y_i)/m) and stipulating that its argument should
4112 * always be integral in the context tableau, i.e., m q = c + \sum a_i y_i.
4113 * Since q is expressed in the tableau as
4114 * c + \sum a_i y_i - m q >= 0
4115 * -c - \sum a_i y_i + m q + m - 1 >= 0
4116 * it is sufficient to add the inequality
4117 * -c - \sum a_i y_i + m q >= 0
4118 * In the part of the context where this inequality does not hold, the
4119 * main tableau is marked as being empty.
4121 static void find_solutions(struct isl_sol *sol, struct isl_tab *tab)
4123 struct isl_context *context;
4124 int r;
4126 if (!tab || sol->error)
4127 goto error;
4129 context = sol->context;
4131 if (tab->empty)
4132 goto done;
4133 if (context->op->is_empty(context))
4134 goto done;
4136 for (r = 0; r >= 0 && tab && !tab->empty; r = restore_lexmin(tab)) {
4137 int flags;
4138 int row;
4139 enum isl_tab_row_sign sgn;
4140 int split = -1;
4141 int n_split = 0;
4143 for (row = tab->n_redundant; row < tab->n_row; ++row) {
4144 if (!isl_tab_var_from_row(tab, row)->is_nonneg)
4145 continue;
4146 sgn = row_sign(tab, sol, row);
4147 if (!sgn)
4148 goto error;
4149 tab->row_sign[row] = sgn;
4150 if (sgn == isl_tab_row_any)
4151 n_split++;
4152 if (sgn == isl_tab_row_any && split == -1)
4153 split = row;
4154 if (sgn == isl_tab_row_neg)
4155 break;
4157 if (row < tab->n_row)
4158 continue;
4159 if (split != -1) {
4160 struct isl_vec *ineq;
4161 if (n_split != 1)
4162 split = context->op->best_split(context, tab);
4163 if (split < 0)
4164 goto error;
4165 ineq = get_row_parameter_ineq(tab, split);
4166 if (!ineq)
4167 goto error;
4168 is_strict(ineq);
4169 reset_any_to_unknown(tab);
4170 tab->row_sign[split] = isl_tab_row_pos;
4171 sol_inc_level(sol);
4172 find_in_pos(sol, tab, ineq->el);
4173 tab->row_sign[split] = isl_tab_row_neg;
4174 isl_seq_neg(ineq->el, ineq->el, ineq->size);
4175 isl_int_sub_ui(ineq->el[0], ineq->el[0], 1);
4176 sol_context_add_ineq(sol, ineq->el, 0, 1);
4177 isl_vec_free(ineq);
4178 if (sol->error)
4179 goto error;
4180 continue;
4182 if (tab->rational)
4183 break;
4184 row = first_non_integer_row(tab, &flags);
4185 if (row < 0)
4186 break;
4187 if (ISL_FL_ISSET(flags, I_PAR)) {
4188 if (ISL_FL_ISSET(flags, I_VAR)) {
4189 if (isl_tab_mark_empty(tab) < 0)
4190 goto error;
4191 break;
4193 row = add_cut(tab, row);
4194 } else if (ISL_FL_ISSET(flags, I_VAR)) {
4195 struct isl_vec *div;
4196 struct isl_vec *ineq;
4197 int d;
4198 div = get_row_split_div(tab, row);
4199 if (!div)
4200 goto error;
4201 d = context->op->get_div(context, tab, div);
4202 isl_vec_free(div);
4203 if (d < 0)
4204 goto error;
4205 ineq = ineq_for_div(context->op->peek_basic_set(context), d);
4206 if (!ineq)
4207 goto error;
4208 sol_inc_level(sol);
4209 no_sol_in_strict(sol, tab, ineq);
4210 isl_seq_neg(ineq->el, ineq->el, ineq->size);
4211 sol_context_add_ineq(sol, ineq->el, 1, 1);
4212 isl_vec_free(ineq);
4213 if (sol->error || !context->op->is_ok(context))
4214 goto error;
4215 tab = set_row_cst_to_div(tab, row, d);
4216 if (context->op->is_empty(context))
4217 break;
4218 } else
4219 row = add_parametric_cut(tab, row, context);
4220 if (row < 0)
4221 goto error;
4223 if (r < 0)
4224 goto error;
4225 done:
4226 sol_add(sol, tab);
4227 isl_tab_free(tab);
4228 return;
4229 error:
4230 isl_tab_free(tab);
4231 sol->error = 1;
4234 /* Does "sol" contain a pair of partial solutions that could potentially
4235 * be merged?
4237 * We currently only check that "sol" is not in an error state
4238 * and that there are at least two partial solutions of which the final two
4239 * are defined at the same level.
4241 static int sol_has_mergeable_solutions(struct isl_sol *sol)
4243 if (sol->error)
4244 return 0;
4245 if (!sol->partial)
4246 return 0;
4247 if (!sol->partial->next)
4248 return 0;
4249 return sol->partial->level == sol->partial->next->level;
4252 /* Compute the lexicographic minimum of the set represented by the main
4253 * tableau "tab" within the context "sol->context_tab".
4255 * As a preprocessing step, we first transfer all the purely parametric
4256 * equalities from the main tableau to the context tableau, i.e.,
4257 * parameters that have been pivoted to a row.
4258 * These equalities are ignored by the main algorithm, because the
4259 * corresponding rows may not be marked as being non-negative.
4260 * In parts of the context where the added equality does not hold,
4261 * the main tableau is marked as being empty.
4263 * Before we embark on the actual computation, we save a copy
4264 * of the context. When we return, we check if there are any
4265 * partial solutions that can potentially be merged. If so,
4266 * we perform a rollback to the initial state of the context.
4267 * The merging of partial solutions happens inside calls to
4268 * sol_dec_level that are pushed onto the undo stack of the context.
4269 * If there are no partial solutions that can potentially be merged
4270 * then the rollback is skipped as it would just be wasted effort.
4272 static void find_solutions_main(struct isl_sol *sol, struct isl_tab *tab)
4274 int row;
4275 void *saved;
4277 if (!tab)
4278 goto error;
4280 sol->level = 0;
4282 for (row = tab->n_redundant; row < tab->n_row; ++row) {
4283 int p;
4284 struct isl_vec *eq;
4286 if (!row_is_parameter_var(tab, row))
4287 continue;
4288 if (tab->row_var[row] < tab->n_param)
4289 p = tab->row_var[row];
4290 else
4291 p = tab->row_var[row]
4292 + tab->n_param - (tab->n_var - tab->n_div);
4294 eq = isl_vec_alloc(tab->mat->ctx, 1+tab->n_param+tab->n_div);
4295 if (!eq)
4296 goto error;
4297 get_row_parameter_line(tab, row, eq->el);
4298 isl_int_neg(eq->el[1 + p], tab->mat->row[row][0]);
4299 eq = isl_vec_normalize(eq);
4301 sol_inc_level(sol);
4302 no_sol_in_strict(sol, tab, eq);
4304 isl_seq_neg(eq->el, eq->el, eq->size);
4305 sol_inc_level(sol);
4306 no_sol_in_strict(sol, tab, eq);
4307 isl_seq_neg(eq->el, eq->el, eq->size);
4309 sol_context_add_eq(sol, eq->el, 1, 1);
4311 isl_vec_free(eq);
4313 if (isl_tab_mark_redundant(tab, row) < 0)
4314 goto error;
4316 if (sol->context->op->is_empty(sol->context))
4317 break;
4319 row = tab->n_redundant - 1;
4322 saved = sol->context->op->save(sol->context);
4324 find_solutions(sol, tab);
4326 if (sol_has_mergeable_solutions(sol))
4327 sol->context->op->restore(sol->context, saved);
4328 else
4329 sol->context->op->discard(saved);
4331 sol->level = 0;
4332 sol_pop(sol);
4334 return;
4335 error:
4336 isl_tab_free(tab);
4337 sol->error = 1;
4340 /* Is the local variable "div" of "bmap" an integer division
4341 * with a known expression that does not involve the "n" variables
4342 * starting at "first"?
4344 static isl_bool is_known_div_not_involving(__isl_keep isl_basic_map *bmap,
4345 unsigned div, unsigned first, unsigned n)
4347 isl_bool unknown, involves;
4349 unknown = isl_basic_map_div_is_marked_unknown(bmap, div);
4350 if (unknown < 0 || unknown)
4351 return isl_bool_not(unknown);
4352 involves = isl_basic_map_div_expr_involves_vars(bmap, div, first, n);
4353 return isl_bool_not(involves);
4356 /* Check if integer division "div" of "src" also occurs in "dst",
4357 * where the integer division "div" is known to involve
4358 * only the first "n_shared" variables.
4359 * If so, return its position within the local variables.
4360 * Otherwise, return a position beyond the local variables.
4362 static isl_size find_div_involving_only(__isl_keep isl_basic_map *dst,
4363 __isl_keep isl_basic_map *src, unsigned div, unsigned n_shared)
4365 int i;
4366 isl_size total;
4367 isl_size n_div;
4369 total = isl_basic_map_dim(dst, isl_dim_all);
4370 n_div = isl_basic_map_dim(dst, isl_dim_div);
4371 if (total < 0 || n_div < 0)
4372 return isl_size_error;
4374 for (i = 0; i < n_div; ++i) {
4375 isl_bool ok;
4377 ok = is_known_div_not_involving(dst, i,
4378 n_shared, total - n_shared);
4379 if (ok < 0)
4380 return isl_size_error;
4381 if (!ok)
4382 continue;
4383 if (isl_seq_eq(dst->div[i], src->div[div], 2 + n_shared))
4384 return i;
4386 return n_div;
4389 /* Check if integer division "div" of "dom" also occurs in "bmap".
4390 * If so, return its position within the divs.
4391 * Otherwise, return a position beyond the integer divisions.
4393 static isl_size find_context_div(__isl_keep isl_basic_map *bmap,
4394 __isl_keep isl_basic_set *dom, unsigned div)
4396 isl_size d_v_div;
4397 isl_size n_div, d_n_div;
4398 isl_bool ok;
4400 d_v_div = isl_basic_set_var_offset(dom, isl_dim_div);
4401 n_div = isl_basic_map_dim(bmap, isl_dim_div);
4402 d_n_div = isl_basic_set_dim(dom, isl_dim_div);
4403 if (d_v_div < 0 || n_div < 0 || d_n_div < 0)
4404 return isl_size_error;
4406 ok = is_known_div_not_involving(bset_to_bmap(dom), div,
4407 d_v_div, d_n_div);
4408 if (ok < 0)
4409 return isl_size_error;
4410 if (!ok)
4411 return n_div;
4413 return find_div_involving_only(bmap, bset_to_bmap(dom), div, d_v_div);
4416 /* Copy integer division "div" of "bmap", which is known to only involve
4417 * the first "n_shared" variables, to "dom" at position "dom_div".
4419 static __isl_give isl_basic_set *copy_div(__isl_take isl_basic_set *dom,
4420 __isl_keep isl_basic_map *bmap,
4421 unsigned div, unsigned n_shared, unsigned dom_div)
4423 isl_vec *v;
4424 isl_size total;
4426 total = isl_basic_set_dim(dom, isl_dim_all);
4427 if (total < 0)
4428 return isl_basic_set_free(dom);
4430 v = isl_vec_alloc(isl_basic_set_get_ctx(dom), 1 + 1 + total);
4431 if (!v)
4432 return isl_basic_set_free(dom);
4434 isl_seq_cpy(v->el, bmap->div[div], 1 + 1 + n_shared);
4435 isl_seq_clr(v->el + 1 + 1 + n_shared, total - n_shared);
4436 dom = isl_basic_set_insert_div(dom, dom_div, v);
4437 dom = isl_basic_set_add_div_constraints(dom, dom_div);
4438 isl_vec_free(v);
4440 return dom;
4443 /* Copy the integer divisions of "bmap" that only involve variables in "dom" and
4444 * that do not already appear in "dom" to "dom".
4446 static __isl_give isl_basic_set *copy_divs(__isl_take isl_basic_set *dom,
4447 __isl_keep isl_basic_map *bmap)
4449 int i;
4450 isl_size dom_n_div, bmap_n_div, total;
4451 isl_size v_out;
4453 dom_n_div = isl_basic_set_dim(dom, isl_dim_div);
4454 bmap_n_div = isl_basic_map_dim(bmap, isl_dim_div);
4455 total = isl_basic_map_dim(bmap, isl_dim_all);
4456 v_out = isl_basic_map_var_offset(bmap, isl_dim_out);
4457 if (dom_n_div < 0 || bmap_n_div < 0 || total < 0 || v_out < 0)
4458 return isl_basic_set_free(dom);
4460 for (i = 0; i < bmap_n_div; ++i) {
4461 isl_bool ok;
4462 isl_size pos;
4464 ok = is_known_div_not_involving(bmap, i, v_out, total - v_out);
4465 if (ok < 0)
4466 return isl_basic_set_free(dom);
4467 if (!ok)
4468 continue;
4469 pos = find_div_involving_only(bset_to_bmap(dom),
4470 bmap, i, v_out);
4471 if (pos < 0)
4472 return isl_basic_set_free(dom);
4473 if (pos < dom_n_div)
4474 continue;
4475 dom = copy_div(dom, bmap, i, v_out, dom_n_div++);
4478 return dom;
4481 /* The correspondence between the variables in the main tableau,
4482 * the context tableau, and the input map and domain is as follows.
4483 * The first n_param and the last n_div variables of the main tableau
4484 * form the variables of the context tableau.
4485 * In the basic map, these n_param variables correspond to the
4486 * parameters and the input dimensions. In the domain, they correspond
4487 * to the parameters and the set dimensions.
4488 * The n_div variables correspond to the integer divisions in the domain.
4489 * To ensure that everything lines up, we may need to copy some of the
4490 * integer divisions of the domain to the map. These have to be placed
4491 * in the same order as those in the context and they have to be placed
4492 * after any other integer divisions that the map may have.
4493 * This function performs the required reordering.
4495 static __isl_give isl_basic_map *align_context_divs(
4496 __isl_take isl_basic_map *bmap, __isl_keep isl_basic_set *dom)
4498 int i;
4499 int common = 0;
4500 int other;
4501 isl_size bmap_n_div, dom_n_div;
4503 bmap_n_div = isl_basic_map_dim(bmap, isl_dim_div);
4504 dom_n_div = isl_basic_set_dim(dom, isl_dim_div);
4505 if (bmap_n_div < 0 || dom_n_div < 0)
4506 return isl_basic_map_free(bmap);
4508 for (i = 0; i < dom_n_div; ++i) {
4509 isl_size pos;
4511 pos = find_context_div(bmap, dom, i);
4512 if (pos < 0)
4513 return isl_basic_map_free(bmap);
4514 if (pos < bmap_n_div)
4515 common++;
4517 other = bmap_n_div - common;
4518 if (dom_n_div - common > 0) {
4519 bmap = isl_basic_map_cow(bmap);
4520 bmap = isl_basic_map_extend(bmap, dom_n_div - common, 0, 0);
4521 if (!bmap)
4522 return NULL;
4524 for (i = 0; i < dom_n_div; ++i) {
4525 isl_size pos = find_context_div(bmap, dom, i);
4526 if (pos < 0)
4527 bmap = isl_basic_map_free(bmap);
4528 if (pos >= bmap_n_div) {
4529 pos = isl_basic_map_alloc_div(bmap);
4530 if (pos < 0)
4531 goto error;
4532 isl_int_set_si(bmap->div[pos][0], 0);
4533 bmap_n_div++;
4535 if (pos != other + i)
4536 bmap = isl_basic_map_swap_div(bmap, pos, other + i);
4538 return bmap;
4539 error:
4540 isl_basic_map_free(bmap);
4541 return NULL;
4544 /* Base case of isl_tab_basic_map_partial_lexopt, after removing
4545 * some obvious symmetries.
4547 * We make sure the divs in the domain are properly ordered,
4548 * because they will be added one by one in the given order
4549 * during the construction of the solution map.
4550 * Furthermore, make sure that the known integer divisions
4551 * appear before any unknown integer division because the solution
4552 * may depend on the known integer divisions, while anything that
4553 * depends on any variable starting from the first unknown integer
4554 * division is ignored in sol_pma_add.
4555 * First copy over any integer divisions from "bmap" that do not
4556 * already appear in "dom". This ensures that the tableau
4557 * will not be split on the corresponding integer division constraints
4558 * since they will be known to hold in "dom".
4560 static struct isl_sol *basic_map_partial_lexopt_base_sol(
4561 __isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom,
4562 __isl_give isl_set **empty, int max,
4563 struct isl_sol *(*init)(__isl_keep isl_basic_map *bmap,
4564 __isl_take isl_basic_set *dom, int track_empty, int max))
4566 struct isl_tab *tab;
4567 struct isl_sol *sol = NULL;
4568 struct isl_context *context;
4570 dom = copy_divs(dom, bmap);
4571 dom = isl_basic_set_sort_divs(dom);
4572 bmap = align_context_divs(bmap, dom);
4573 sol = init(bmap, dom, !!empty, max);
4574 if (!sol)
4575 goto error;
4577 context = sol->context;
4578 if (isl_basic_set_plain_is_empty(context->op->peek_basic_set(context)))
4579 /* nothing */;
4580 else if (isl_basic_map_plain_is_empty(bmap)) {
4581 if (sol->add_empty)
4582 sol->add_empty(sol,
4583 isl_basic_set_copy(context->op->peek_basic_set(context)));
4584 } else {
4585 tab = tab_for_lexmin(bmap,
4586 context->op->peek_basic_set(context), 1, max);
4587 tab = context->op->detect_nonnegative_parameters(context, tab);
4588 find_solutions_main(sol, tab);
4590 if (sol->error)
4591 goto error;
4593 isl_basic_map_free(bmap);
4594 return sol;
4595 error:
4596 sol_free(sol);
4597 isl_basic_map_free(bmap);
4598 return NULL;
4601 /* Base case of isl_tab_basic_map_partial_lexopt, after removing
4602 * some obvious symmetries.
4604 * We call basic_map_partial_lexopt_base_sol and extract the results.
4606 static __isl_give isl_map *basic_map_partial_lexopt_base(
4607 __isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom,
4608 __isl_give isl_set **empty, int max)
4610 isl_map *result = NULL;
4611 struct isl_sol *sol;
4612 struct isl_sol_map *sol_map;
4614 sol = basic_map_partial_lexopt_base_sol(bmap, dom, empty, max,
4615 &sol_map_init);
4616 if (!sol)
4617 return NULL;
4618 sol_map = (struct isl_sol_map *) sol;
4620 result = isl_map_copy(sol_map->map);
4621 if (empty)
4622 *empty = isl_set_copy(sol_map->empty);
4623 sol_free(&sol_map->sol);
4624 return result;
4627 /* Return a count of the number of occurrences of the "n" first
4628 * variables in the inequality constraints of "bmap".
4630 static __isl_give int *count_occurrences(__isl_keep isl_basic_map *bmap,
4631 int n)
4633 int i, j;
4634 isl_ctx *ctx;
4635 int *occurrences;
4637 if (!bmap)
4638 return NULL;
4639 ctx = isl_basic_map_get_ctx(bmap);
4640 occurrences = isl_calloc_array(ctx, int, n);
4641 if (!occurrences)
4642 return NULL;
4644 for (i = 0; i < bmap->n_ineq; ++i) {
4645 for (j = 0; j < n; ++j) {
4646 if (!isl_int_is_zero(bmap->ineq[i][1 + j]))
4647 occurrences[j]++;
4651 return occurrences;
4654 /* Do all of the "n" variables with non-zero coefficients in "c"
4655 * occur in exactly a single constraint.
4656 * "occurrences" is an array of length "n" containing the number
4657 * of occurrences of each of the variables in the inequality constraints.
4659 static int single_occurrence(int n, isl_int *c, int *occurrences)
4661 int i;
4663 for (i = 0; i < n; ++i) {
4664 if (isl_int_is_zero(c[i]))
4665 continue;
4666 if (occurrences[i] != 1)
4667 return 0;
4670 return 1;
4673 /* Do all of the "n" initial variables that occur in inequality constraint
4674 * "ineq" of "bmap" only occur in that constraint?
4676 static int all_single_occurrence(__isl_keep isl_basic_map *bmap, int ineq,
4677 int n)
4679 int i, j;
4681 for (i = 0; i < n; ++i) {
4682 if (isl_int_is_zero(bmap->ineq[ineq][1 + i]))
4683 continue;
4684 for (j = 0; j < bmap->n_ineq; ++j) {
4685 if (j == ineq)
4686 continue;
4687 if (!isl_int_is_zero(bmap->ineq[j][1 + i]))
4688 return 0;
4692 return 1;
4695 /* Structure used during detection of parallel constraints.
4696 * n_in: number of "input" variables: isl_dim_param + isl_dim_in
4697 * n_out: number of "output" variables: isl_dim_out + isl_dim_div
4698 * val: the coefficients of the output variables
4700 struct isl_constraint_equal_info {
4701 unsigned n_in;
4702 unsigned n_out;
4703 isl_int *val;
4706 /* Check whether the coefficients of the output variables
4707 * of the constraint in "entry" are equal to info->val.
4709 static isl_bool constraint_equal(const void *entry, const void *val)
4711 isl_int **row = (isl_int **)entry;
4712 const struct isl_constraint_equal_info *info = val;
4713 int eq;
4715 eq = isl_seq_eq((*row) + 1 + info->n_in, info->val, info->n_out);
4716 return isl_bool_ok(eq);
4719 /* Check whether "bmap" has a pair of constraints that have
4720 * the same coefficients for the output variables.
4721 * Note that the coefficients of the existentially quantified
4722 * variables need to be zero since the existentially quantified
4723 * of the result are usually not the same as those of the input.
4724 * Furthermore, check that each of the input variables that occur
4725 * in those constraints does not occur in any other constraint.
4726 * If so, return true and return the row indices of the two constraints
4727 * in *first and *second.
4729 static isl_bool parallel_constraints(__isl_keep isl_basic_map *bmap,
4730 int *first, int *second)
4732 int i;
4733 isl_ctx *ctx;
4734 int *occurrences = NULL;
4735 struct isl_hash_table *table = NULL;
4736 struct isl_hash_table_entry *entry;
4737 struct isl_constraint_equal_info info;
4738 isl_size nparam, n_in, n_out, n_div;
4740 ctx = isl_basic_map_get_ctx(bmap);
4741 table = isl_hash_table_alloc(ctx, bmap->n_ineq);
4742 if (!table)
4743 goto error;
4745 nparam = isl_basic_map_dim(bmap, isl_dim_param);
4746 n_in = isl_basic_map_dim(bmap, isl_dim_in);
4747 n_out = isl_basic_map_dim(bmap, isl_dim_out);
4748 n_div = isl_basic_map_dim(bmap, isl_dim_div);
4749 if (nparam < 0 || n_in < 0 || n_out < 0 || n_div < 0)
4750 goto error;
4751 info.n_in = nparam + n_in;
4752 occurrences = count_occurrences(bmap, info.n_in);
4753 if (info.n_in && !occurrences)
4754 goto error;
4755 info.n_out = n_out + n_div;
4756 for (i = 0; i < bmap->n_ineq; ++i) {
4757 uint32_t hash;
4759 info.val = bmap->ineq[i] + 1 + info.n_in;
4760 if (!isl_seq_any_non_zero(info.val, n_out))
4761 continue;
4762 if (isl_seq_any_non_zero(info.val + n_out, n_div))
4763 continue;
4764 if (!single_occurrence(info.n_in, bmap->ineq[i] + 1,
4765 occurrences))
4766 continue;
4767 hash = isl_seq_get_hash(info.val, info.n_out);
4768 entry = isl_hash_table_find(ctx, table, hash,
4769 constraint_equal, &info, 1);
4770 if (!entry)
4771 goto error;
4772 if (entry->data)
4773 break;
4774 entry->data = &bmap->ineq[i];
4777 if (i < bmap->n_ineq) {
4778 *first = ((isl_int **)entry->data) - bmap->ineq;
4779 *second = i;
4782 isl_hash_table_free(ctx, table);
4783 free(occurrences);
4785 return isl_bool_ok(i < bmap->n_ineq);
4786 error:
4787 isl_hash_table_free(ctx, table);
4788 free(occurrences);
4789 return isl_bool_error;
4792 /* Given a set of upper bounds in "var", add constraints to "bset"
4793 * that make the i-th bound smallest.
4795 * In particular, if there are n bounds b_i, then add the constraints
4797 * b_i <= b_j for j > i
4798 * b_i < b_j for j < i
4800 static __isl_give isl_basic_set *select_minimum(__isl_take isl_basic_set *bset,
4801 __isl_keep isl_mat *var, int i)
4803 isl_ctx *ctx;
4804 int j, k;
4806 ctx = isl_mat_get_ctx(var);
4808 for (j = 0; j < var->n_row; ++j) {
4809 if (j == i)
4810 continue;
4811 k = isl_basic_set_alloc_inequality(bset);
4812 if (k < 0)
4813 goto error;
4814 isl_seq_combine(bset->ineq[k], ctx->one, var->row[j],
4815 ctx->negone, var->row[i], var->n_col);
4816 isl_int_set_si(bset->ineq[k][var->n_col], 0);
4817 if (j < i)
4818 isl_int_sub_ui(bset->ineq[k][0], bset->ineq[k][0], 1);
4821 bset = isl_basic_set_finalize(bset);
4823 return bset;
4824 error:
4825 isl_basic_set_free(bset);
4826 return NULL;
4829 /* Given a set of upper bounds on the last "input" variable m,
4830 * construct a set that assigns the minimal upper bound to m, i.e.,
4831 * construct a set that divides the space into cells where one
4832 * of the upper bounds is smaller than all the others and assign
4833 * this upper bound to m.
4835 * In particular, if there are n bounds b_i, then the result
4836 * consists of n basic sets, each one of the form
4838 * m = b_i
4839 * b_i <= b_j for j > i
4840 * b_i < b_j for j < i
4842 static __isl_give isl_set *set_minimum(__isl_take isl_space *space,
4843 __isl_take isl_mat *var)
4845 int i, k;
4846 isl_basic_set *bset = NULL;
4847 isl_set *set = NULL;
4849 if (!space || !var)
4850 goto error;
4852 set = isl_set_alloc_space(isl_space_copy(space),
4853 var->n_row, ISL_SET_DISJOINT);
4855 for (i = 0; i < var->n_row; ++i) {
4856 bset = isl_basic_set_alloc_space(isl_space_copy(space), 0,
4857 1, var->n_row - 1);
4858 k = isl_basic_set_alloc_equality(bset);
4859 if (k < 0)
4860 goto error;
4861 isl_seq_cpy(bset->eq[k], var->row[i], var->n_col);
4862 isl_int_set_si(bset->eq[k][var->n_col], -1);
4863 bset = select_minimum(bset, var, i);
4864 set = isl_set_add_basic_set(set, bset);
4867 isl_space_free(space);
4868 isl_mat_free(var);
4869 return set;
4870 error:
4871 isl_basic_set_free(bset);
4872 isl_set_free(set);
4873 isl_space_free(space);
4874 isl_mat_free(var);
4875 return NULL;
4878 /* Given that the last input variable of "bmap" represents the minimum
4879 * of the bounds in "cst", check whether we need to split the domain
4880 * based on which bound attains the minimum.
4882 * A split is needed when the minimum appears in an integer division
4883 * or in an equality. Otherwise, it is only needed if it appears in
4884 * an upper bound that is different from the upper bounds on which it
4885 * is defined.
4887 static isl_bool need_split_basic_map(__isl_keep isl_basic_map *bmap,
4888 __isl_keep isl_mat *cst)
4890 int i, j;
4891 isl_bool involves;
4892 isl_size total;
4893 unsigned pos;
4895 pos = cst->n_col - 1;
4896 total = isl_basic_map_dim(bmap, isl_dim_all);
4897 if (total < 0)
4898 return isl_bool_error;
4900 involves = isl_basic_map_any_div_involves_vars(bmap, pos, 1);
4901 if (involves < 0 || involves)
4902 return involves;
4904 for (i = 0; i < bmap->n_eq; ++i)
4905 if (!isl_int_is_zero(bmap->eq[i][1 + pos]))
4906 return isl_bool_true;
4908 for (i = 0; i < bmap->n_ineq; ++i) {
4909 if (isl_int_is_nonneg(bmap->ineq[i][1 + pos]))
4910 continue;
4911 if (!isl_int_is_negone(bmap->ineq[i][1 + pos]))
4912 return isl_bool_true;
4913 if (isl_seq_any_non_zero(bmap->ineq[i] + 1 + pos + 1,
4914 total - pos - 1))
4915 return isl_bool_true;
4917 for (j = 0; j < cst->n_row; ++j)
4918 if (isl_seq_eq(bmap->ineq[i], cst->row[j], cst->n_col))
4919 break;
4920 if (j >= cst->n_row)
4921 return isl_bool_true;
4924 return isl_bool_false;
4927 /* Given that the last set variable of "bset" represents the minimum
4928 * of the bounds in "cst", check whether we need to split the domain
4929 * based on which bound attains the minimum.
4931 * We simply call need_split_basic_map here. This is safe because
4932 * the position of the minimum is computed from "cst" and not
4933 * from "bmap".
4935 static isl_bool need_split_basic_set(__isl_keep isl_basic_set *bset,
4936 __isl_keep isl_mat *cst)
4938 return need_split_basic_map(bset_to_bmap(bset), cst);
4941 /* Given that the last set variable of "set" represents the minimum
4942 * of the bounds in "cst", check whether we need to split the domain
4943 * based on which bound attains the minimum.
4945 static isl_bool need_split_set(__isl_keep isl_set *set, __isl_keep isl_mat *cst)
4947 int i;
4949 for (i = 0; i < set->n; ++i) {
4950 isl_bool split;
4952 split = need_split_basic_set(set->p[i], cst);
4953 if (split < 0 || split)
4954 return split;
4957 return isl_bool_false;
4960 /* Given a map of which the last input variable is the minimum
4961 * of the bounds in "cst", split each basic set in the set
4962 * in pieces where one of the bounds is (strictly) smaller than the others.
4963 * This subdivision is given in "min_expr".
4964 * The variable is subsequently projected out.
4966 * We only do the split when it is needed.
4967 * For example if the last input variable m = min(a,b) and the only
4968 * constraints in the given basic set are lower bounds on m,
4969 * i.e., l <= m = min(a,b), then we can simply project out m
4970 * to obtain l <= a and l <= b, without having to split on whether
4971 * m is equal to a or b.
4973 static __isl_give isl_map *split_domain(__isl_take isl_map *opt,
4974 __isl_take isl_set *min_expr, __isl_take isl_mat *cst)
4976 isl_size n_in;
4977 int i;
4978 isl_space *space;
4979 isl_map *res;
4981 n_in = isl_map_dim(opt, isl_dim_in);
4982 if (n_in < 0 || !min_expr || !cst)
4983 goto error;
4985 space = isl_map_get_space(opt);
4986 space = isl_space_drop_dims(space, isl_dim_in, n_in - 1, 1);
4987 res = isl_map_empty(space);
4989 for (i = 0; i < opt->n; ++i) {
4990 isl_map *map;
4991 isl_bool split;
4993 map = isl_map_from_basic_map(isl_basic_map_copy(opt->p[i]));
4994 split = need_split_basic_map(opt->p[i], cst);
4995 if (split < 0)
4996 map = isl_map_free(map);
4997 else if (split)
4998 map = isl_map_intersect_domain(map,
4999 isl_set_copy(min_expr));
5000 map = isl_map_remove_dims(map, isl_dim_in, n_in - 1, 1);
5002 res = isl_map_union_disjoint(res, map);
5005 isl_map_free(opt);
5006 isl_set_free(min_expr);
5007 isl_mat_free(cst);
5008 return res;
5009 error:
5010 isl_map_free(opt);
5011 isl_set_free(min_expr);
5012 isl_mat_free(cst);
5013 return NULL;
5016 /* Given a set of which the last set variable is the minimum
5017 * of the bounds in "cst", split each basic set in the set
5018 * in pieces where one of the bounds is (strictly) smaller than the others.
5019 * This subdivision is given in "min_expr".
5020 * The variable is subsequently projected out.
5022 static __isl_give isl_set *split(__isl_take isl_set *empty,
5023 __isl_take isl_set *min_expr, __isl_take isl_mat *cst)
5025 isl_map *map;
5027 map = isl_map_from_domain(empty);
5028 map = split_domain(map, min_expr, cst);
5029 empty = isl_map_domain(map);
5031 return empty;
5034 static __isl_give isl_map *basic_map_partial_lexopt(
5035 __isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom,
5036 __isl_give isl_set **empty, int max);
5038 /* This function is called from basic_map_partial_lexopt_symm.
5039 * The last variable of "bmap" and "dom" corresponds to the minimum
5040 * of the bounds in "cst". "map_space" is the space of the original
5041 * input relation (of basic_map_partial_lexopt_symm) and "set_space"
5042 * is the space of the original domain.
5044 * We recursively call basic_map_partial_lexopt and then plug in
5045 * the definition of the minimum in the result.
5047 static __isl_give isl_map *basic_map_partial_lexopt_symm_core(
5048 __isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom,
5049 __isl_give isl_set **empty, int max, __isl_take isl_mat *cst,
5050 __isl_take isl_space *map_space, __isl_take isl_space *set_space)
5052 isl_map *opt;
5053 isl_set *min_expr;
5055 min_expr = set_minimum(isl_basic_set_get_space(dom), isl_mat_copy(cst));
5057 opt = basic_map_partial_lexopt(bmap, dom, empty, max);
5059 if (empty) {
5060 *empty = split(*empty,
5061 isl_set_copy(min_expr), isl_mat_copy(cst));
5062 *empty = isl_set_reset_space(*empty, set_space);
5065 opt = split_domain(opt, min_expr, cst);
5066 opt = isl_map_reset_space(opt, map_space);
5068 return opt;
5071 /* Extract a domain from "bmap" for the purpose of computing
5072 * a lexicographic optimum.
5074 * This function is only called when the caller wants to compute a full
5075 * lexicographic optimum, i.e., without specifying a domain. In this case,
5076 * the caller is not interested in the part of the domain space where
5077 * there is no solution and the domain can be initialized to those constraints
5078 * of "bmap" that only involve the parameters and the input dimensions.
5079 * This relieves the parametric programming engine from detecting those
5080 * inequalities and transferring them to the context. More importantly,
5081 * it ensures that those inequalities are transferred first and not
5082 * intermixed with inequalities that actually split the domain.
5084 * If the caller does not require the absence of existentially quantified
5085 * variables in the result (i.e., if ISL_OPT_QE is not set in "flags"),
5086 * then the actual domain of "bmap" can be used. This ensures that
5087 * the domain does not need to be split at all just to separate out
5088 * pieces of the domain that do not have a solution from piece that do.
5089 * This domain cannot be used in general because it may involve
5090 * (unknown) existentially quantified variables which will then also
5091 * appear in the solution.
5093 static __isl_give isl_basic_set *extract_domain(__isl_keep isl_basic_map *bmap,
5094 unsigned flags)
5096 isl_size n_div;
5097 isl_size n_out;
5099 n_div = isl_basic_map_dim(bmap, isl_dim_div);
5100 n_out = isl_basic_map_dim(bmap, isl_dim_out);
5101 if (n_div < 0 || n_out < 0)
5102 return NULL;
5103 bmap = isl_basic_map_copy(bmap);
5104 if (ISL_FL_ISSET(flags, ISL_OPT_QE)) {
5105 bmap = isl_basic_map_drop_constraints_involving_dims(bmap,
5106 isl_dim_div, 0, n_div);
5107 bmap = isl_basic_map_drop_constraints_involving_dims(bmap,
5108 isl_dim_out, 0, n_out);
5110 return isl_basic_map_domain(bmap);
5113 #undef TYPE
5114 #define TYPE isl_map
5115 #undef SUFFIX
5116 #define SUFFIX
5117 #define isl_map_pullback_multi_aff isl_map_preimage_domain_multi_aff
5118 #include "isl_tab_lexopt_templ.c"
5120 /* Extract the subsequence of the sample value of "tab"
5121 * starting at "pos" and of length "len".
5123 static __isl_give isl_vec *extract_sample_sequence(struct isl_tab *tab,
5124 int pos, int len)
5126 int i;
5127 isl_ctx *ctx;
5128 isl_vec *v;
5130 ctx = isl_tab_get_ctx(tab);
5131 v = isl_vec_alloc(ctx, len);
5132 if (!v)
5133 return NULL;
5134 for (i = 0; i < len; ++i) {
5135 if (!tab->var[pos + i].is_row) {
5136 isl_int_set_si(v->el[i], 0);
5137 } else {
5138 int row;
5140 row = tab->var[pos + i].index;
5141 isl_int_divexact(v->el[i], tab->mat->row[row][1],
5142 tab->mat->row[row][0]);
5146 return v;
5149 /* Check if the sequence of variables starting at "pos"
5150 * represents a trivial solution according to "trivial".
5151 * That is, is the result of applying "trivial" to this sequence
5152 * equal to the zero vector?
5154 static isl_bool region_is_trivial(struct isl_tab *tab, int pos,
5155 __isl_keep isl_mat *trivial)
5157 isl_size n, len;
5158 isl_vec *v;
5159 isl_bool is_trivial;
5161 n = isl_mat_rows(trivial);
5162 if (n < 0)
5163 return isl_bool_error;
5165 if (n == 0)
5166 return isl_bool_false;
5168 len = isl_mat_cols(trivial);
5169 if (len < 0)
5170 return isl_bool_error;
5171 v = extract_sample_sequence(tab, pos, len);
5172 v = isl_mat_vec_product(isl_mat_copy(trivial), v);
5173 is_trivial = isl_vec_is_zero(v);
5174 isl_vec_free(v);
5176 return is_trivial;
5179 /* Global internal data for isl_tab_basic_set_non_trivial_lexmin.
5181 * "n_op" is the number of initial coordinates to optimize,
5182 * as passed to isl_tab_basic_set_non_trivial_lexmin.
5183 * "region" is the "n_region"-sized array of regions passed
5184 * to isl_tab_basic_set_non_trivial_lexmin.
5186 * "tab" is the tableau that corresponds to the ILP problem.
5187 * "local" is an array of local data structure, one for each
5188 * (potential) level of the backtracking procedure of
5189 * isl_tab_basic_set_non_trivial_lexmin.
5190 * "v" is a pre-allocated vector that can be used for adding
5191 * constraints to the tableau.
5193 * "sol" contains the best solution found so far.
5194 * It is initialized to a vector of size zero.
5196 struct isl_lexmin_data {
5197 int n_op;
5198 int n_region;
5199 struct isl_trivial_region *region;
5201 struct isl_tab *tab;
5202 struct isl_local_region *local;
5203 isl_vec *v;
5205 isl_vec *sol;
5208 /* Return the index of the first trivial region, "n_region" if all regions
5209 * are non-trivial or -1 in case of error.
5211 static int first_trivial_region(struct isl_lexmin_data *data)
5213 int i;
5215 for (i = 0; i < data->n_region; ++i) {
5216 isl_bool trivial;
5217 trivial = region_is_trivial(data->tab, data->region[i].pos,
5218 data->region[i].trivial);
5219 if (trivial < 0)
5220 return -1;
5221 if (trivial)
5222 return i;
5225 return data->n_region;
5228 /* Check if the solution is optimal, i.e., whether the first
5229 * n_op entries are zero.
5231 static int is_optimal(__isl_keep isl_vec *sol, int n_op)
5233 int i;
5235 for (i = 0; i < n_op; ++i)
5236 if (!isl_int_is_zero(sol->el[1 + i]))
5237 return 0;
5238 return 1;
5241 /* Add constraints to "tab" that ensure that any solution is significantly
5242 * better than that represented by "sol". That is, find the first
5243 * relevant (within first n_op) non-zero coefficient and force it (along
5244 * with all previous coefficients) to be zero.
5245 * If the solution is already optimal (all relevant coefficients are zero),
5246 * then just mark the table as empty.
5247 * "n_zero" is the number of coefficients that have been forced zero
5248 * by previous calls to this function at the same level.
5249 * Return the updated number of forced zero coefficients or -1 on error.
5251 * This function assumes that at least 2 * (n_op - n_zero) more rows and
5252 * at least 2 * (n_op - n_zero) more elements in the constraint array
5253 * are available in the tableau.
5255 static int force_better_solution(struct isl_tab *tab,
5256 __isl_keep isl_vec *sol, int n_op, int n_zero)
5258 int i, n;
5259 isl_ctx *ctx;
5260 isl_vec *v = NULL;
5262 if (!sol)
5263 return -1;
5265 for (i = n_zero; i < n_op; ++i)
5266 if (!isl_int_is_zero(sol->el[1 + i]))
5267 break;
5269 if (i == n_op) {
5270 if (isl_tab_mark_empty(tab) < 0)
5271 return -1;
5272 return n_op;
5275 ctx = isl_vec_get_ctx(sol);
5276 v = isl_vec_alloc(ctx, 1 + tab->n_var);
5277 if (!v)
5278 return -1;
5280 n = i + 1;
5281 for (; i >= n_zero; --i) {
5282 v = isl_vec_clr(v);
5283 isl_int_set_si(v->el[1 + i], -1);
5284 if (add_lexmin_eq(tab, v->el) < 0)
5285 goto error;
5288 isl_vec_free(v);
5289 return n;
5290 error:
5291 isl_vec_free(v);
5292 return -1;
5295 /* Fix triviality direction "dir" of the given region to zero.
5297 * This function assumes that at least two more rows and at least
5298 * two more elements in the constraint array are available in the tableau.
5300 static isl_stat fix_zero(struct isl_tab *tab, struct isl_trivial_region *region,
5301 int dir, struct isl_lexmin_data *data)
5303 isl_size len;
5305 data->v = isl_vec_clr(data->v);
5306 if (!data->v)
5307 return isl_stat_error;
5308 len = isl_mat_cols(region->trivial);
5309 if (len < 0)
5310 return isl_stat_error;
5311 isl_seq_cpy(data->v->el + 1 + region->pos, region->trivial->row[dir],
5312 len);
5313 if (add_lexmin_eq(tab, data->v->el) < 0)
5314 return isl_stat_error;
5316 return isl_stat_ok;
5319 /* This function selects case "side" for non-triviality region "region",
5320 * assuming all the equality constraints have been imposed already.
5321 * In particular, the triviality direction side/2 is made positive
5322 * if side is even and made negative if side is odd.
5324 * This function assumes that at least one more row and at least
5325 * one more element in the constraint array are available in the tableau.
5327 static struct isl_tab *pos_neg(struct isl_tab *tab,
5328 struct isl_trivial_region *region,
5329 int side, struct isl_lexmin_data *data)
5331 isl_size len;
5333 data->v = isl_vec_clr(data->v);
5334 if (!data->v)
5335 goto error;
5336 isl_int_set_si(data->v->el[0], -1);
5337 len = isl_mat_cols(region->trivial);
5338 if (len < 0)
5339 goto error;
5340 if (side % 2 == 0)
5341 isl_seq_cpy(data->v->el + 1 + region->pos,
5342 region->trivial->row[side / 2], len);
5343 else
5344 isl_seq_neg(data->v->el + 1 + region->pos,
5345 region->trivial->row[side / 2], len);
5346 return add_lexmin_ineq(tab, data->v->el);
5347 error:
5348 isl_tab_free(tab);
5349 return NULL;
5352 /* Local data at each level of the backtracking procedure of
5353 * isl_tab_basic_set_non_trivial_lexmin.
5355 * "update" is set if a solution has been found in the current case
5356 * of this level, such that a better solution needs to be enforced
5357 * in the next case.
5358 * "n_zero" is the number of initial coordinates that have already
5359 * been forced to be zero at this level.
5360 * "region" is the non-triviality region considered at this level.
5361 * "side" is the index of the current case at this level.
5362 * "n" is the number of triviality directions.
5363 * "snap" is a snapshot of the tableau holding a state that needs
5364 * to be satisfied by all subsequent cases.
5366 struct isl_local_region {
5367 int update;
5368 int n_zero;
5369 int region;
5370 int side;
5371 int n;
5372 struct isl_tab_undo *snap;
5375 /* Initialize the global data structure "data" used while solving
5376 * the ILP problem "bset".
5378 static isl_stat init_lexmin_data(struct isl_lexmin_data *data,
5379 __isl_keep isl_basic_set *bset)
5381 isl_ctx *ctx;
5383 ctx = isl_basic_set_get_ctx(bset);
5385 data->tab = tab_for_lexmin(bset, NULL, 0, 0);
5386 if (!data->tab)
5387 return isl_stat_error;
5389 data->v = isl_vec_alloc(ctx, 1 + data->tab->n_var);
5390 if (!data->v)
5391 return isl_stat_error;
5392 data->local = isl_calloc_array(ctx, struct isl_local_region,
5393 data->n_region);
5394 if (data->n_region && !data->local)
5395 return isl_stat_error;
5397 data->sol = isl_vec_alloc(ctx, 0);
5399 return isl_stat_ok;
5402 /* Mark all outer levels as requiring a better solution
5403 * in the next cases.
5405 static void update_outer_levels(struct isl_lexmin_data *data, int level)
5407 int i;
5409 for (i = 0; i < level; ++i)
5410 data->local[i].update = 1;
5413 /* Initialize "local" to refer to region "region" and
5414 * to initiate processing at this level.
5416 static isl_stat init_local_region(struct isl_local_region *local, int region,
5417 struct isl_lexmin_data *data)
5419 isl_size n = isl_mat_rows(data->region[region].trivial);
5421 if (n < 0)
5422 return isl_stat_error;
5423 local->n = n;
5424 local->region = region;
5425 local->side = 0;
5426 local->update = 0;
5427 local->n_zero = 0;
5429 return isl_stat_ok;
5432 /* What to do next after entering a level of the backtracking procedure.
5434 * error: some error has occurred; abort
5435 * done: an optimal solution has been found; stop search
5436 * backtrack: backtrack to the previous level
5437 * handle: add the constraints for the current level and
5438 * move to the next level
5440 enum isl_next {
5441 isl_next_error = -1,
5442 isl_next_done,
5443 isl_next_backtrack,
5444 isl_next_handle,
5447 /* Have all cases of the current region been considered?
5448 * If there are n directions, then there are 2n cases.
5450 * The constraints in the current tableau are imposed
5451 * in all subsequent cases. This means that if the current
5452 * tableau is empty, then none of those cases should be considered
5453 * anymore and all cases have effectively been considered.
5455 static int finished_all_cases(struct isl_local_region *local,
5456 struct isl_lexmin_data *data)
5458 if (data->tab->empty)
5459 return 1;
5460 return local->side >= 2 * local->n;
5463 /* Enter level "level" of the backtracking search and figure out
5464 * what to do next. "init" is set if the level was entered
5465 * from a higher level and needs to be initialized.
5466 * Otherwise, the level is entered as a result of backtracking and
5467 * the tableau needs to be restored to a position that can
5468 * be used for the next case at this level.
5469 * The snapshot is assumed to have been saved in the previous case,
5470 * before the constraints specific to that case were added.
5472 * In the initialization case, the local region is initialized
5473 * to point to the first violated region.
5474 * If the constraints of all regions are satisfied by the current
5475 * sample of the tableau, then tell the caller to continue looking
5476 * for a better solution or to stop searching if an optimal solution
5477 * has been found.
5479 * If the tableau is empty or if all cases at the current level
5480 * have been considered, then the caller needs to backtrack as well.
5482 static enum isl_next enter_level(int level, int init,
5483 struct isl_lexmin_data *data)
5485 struct isl_local_region *local = &data->local[level];
5487 if (init) {
5488 int r;
5490 data->tab = cut_to_integer_lexmin(data->tab, CUT_ONE);
5491 if (!data->tab)
5492 return isl_next_error;
5493 if (data->tab->empty)
5494 return isl_next_backtrack;
5495 r = first_trivial_region(data);
5496 if (r < 0)
5497 return isl_next_error;
5498 if (r == data->n_region) {
5499 update_outer_levels(data, level);
5500 isl_vec_free(data->sol);
5501 data->sol = isl_tab_get_sample_value(data->tab);
5502 if (!data->sol)
5503 return isl_next_error;
5504 if (is_optimal(data->sol, data->n_op))
5505 return isl_next_done;
5506 return isl_next_backtrack;
5508 if (level >= data->n_region)
5509 isl_die(isl_vec_get_ctx(data->v), isl_error_internal,
5510 "nesting level too deep",
5511 return isl_next_error);
5512 if (init_local_region(local, r, data) < 0)
5513 return isl_next_error;
5514 if (isl_tab_extend_cons(data->tab,
5515 2 * local->n + 2 * data->n_op) < 0)
5516 return isl_next_error;
5517 } else {
5518 if (isl_tab_rollback(data->tab, local->snap) < 0)
5519 return isl_next_error;
5522 if (finished_all_cases(local, data))
5523 return isl_next_backtrack;
5524 return isl_next_handle;
5527 /* If a solution has been found in the previous case at this level
5528 * (marked by local->update being set), then add constraints
5529 * that enforce a better solution in the present and all following cases.
5530 * The constraints only need to be imposed once because they are
5531 * included in the snapshot (taken in pick_side) that will be used in
5532 * subsequent cases.
5534 static isl_stat better_next_side(struct isl_local_region *local,
5535 struct isl_lexmin_data *data)
5537 if (!local->update)
5538 return isl_stat_ok;
5540 local->n_zero = force_better_solution(data->tab,
5541 data->sol, data->n_op, local->n_zero);
5542 if (local->n_zero < 0)
5543 return isl_stat_error;
5545 local->update = 0;
5547 return isl_stat_ok;
5550 /* Add constraints to data->tab that select the current case (local->side)
5551 * at the current level.
5553 * If the linear combinations v should not be zero, then the cases are
5554 * v_0 >= 1
5555 * v_0 <= -1
5556 * v_0 = 0 and v_1 >= 1
5557 * v_0 = 0 and v_1 <= -1
5558 * v_0 = 0 and v_1 = 0 and v_2 >= 1
5559 * v_0 = 0 and v_1 = 0 and v_2 <= -1
5560 * ...
5561 * in this order.
5563 * A snapshot is taken after the equality constraint (if any) has been added
5564 * such that the next case can start off from this position.
5565 * The rollback to this position is performed in enter_level.
5567 static isl_stat pick_side(struct isl_local_region *local,
5568 struct isl_lexmin_data *data)
5570 struct isl_trivial_region *region;
5571 int side, base;
5573 region = &data->region[local->region];
5574 side = local->side;
5575 base = 2 * (side/2);
5577 if (side == base && base >= 2 &&
5578 fix_zero(data->tab, region, base / 2 - 1, data) < 0)
5579 return isl_stat_error;
5581 local->snap = isl_tab_snap(data->tab);
5582 if (isl_tab_push_basis(data->tab) < 0)
5583 return isl_stat_error;
5585 data->tab = pos_neg(data->tab, region, side, data);
5586 if (!data->tab)
5587 return isl_stat_error;
5588 return isl_stat_ok;
5591 /* Free the memory associated to "data".
5593 static void clear_lexmin_data(struct isl_lexmin_data *data)
5595 free(data->local);
5596 isl_vec_free(data->v);
5597 isl_tab_free(data->tab);
5600 /* Return the lexicographically smallest non-trivial solution of the
5601 * given ILP problem.
5603 * All variables are assumed to be non-negative.
5605 * n_op is the number of initial coordinates to optimize.
5606 * That is, once a solution has been found, we will only continue looking
5607 * for solutions that result in significantly better values for those
5608 * initial coordinates. That is, we only continue looking for solutions
5609 * that increase the number of initial zeros in this sequence.
5611 * A solution is non-trivial, if it is non-trivial on each of the
5612 * specified regions. Each region represents a sequence of
5613 * triviality directions on a sequence of variables that starts
5614 * at a given position. A solution is non-trivial on such a region if
5615 * at least one of the triviality directions is non-zero
5616 * on that sequence of variables.
5618 * Whenever a conflict is encountered, all constraints involved are
5619 * reported to the caller through a call to "conflict".
5621 * We perform a simple branch-and-bound backtracking search.
5622 * Each level in the search represents an initially trivial region
5623 * that is forced to be non-trivial.
5624 * At each level we consider 2 * n cases, where n
5625 * is the number of triviality directions.
5626 * In terms of those n directions v_i, we consider the cases
5627 * v_0 >= 1
5628 * v_0 <= -1
5629 * v_0 = 0 and v_1 >= 1
5630 * v_0 = 0 and v_1 <= -1
5631 * v_0 = 0 and v_1 = 0 and v_2 >= 1
5632 * v_0 = 0 and v_1 = 0 and v_2 <= -1
5633 * ...
5634 * in this order.
5636 __isl_give isl_vec *isl_tab_basic_set_non_trivial_lexmin(
5637 __isl_take isl_basic_set *bset, int n_op, int n_region,
5638 struct isl_trivial_region *region,
5639 int (*conflict)(int con, void *user), void *user)
5641 struct isl_lexmin_data data = { n_op, n_region, region };
5642 int level, init;
5644 if (!bset)
5645 return NULL;
5647 if (init_lexmin_data(&data, bset) < 0)
5648 goto error;
5649 data.tab->conflict = conflict;
5650 data.tab->conflict_user = user;
5652 level = 0;
5653 init = 1;
5655 while (level >= 0) {
5656 enum isl_next next;
5657 struct isl_local_region *local = &data.local[level];
5659 next = enter_level(level, init, &data);
5660 if (next < 0)
5661 goto error;
5662 if (next == isl_next_done)
5663 break;
5664 if (next == isl_next_backtrack) {
5665 level--;
5666 init = 0;
5667 continue;
5670 if (better_next_side(local, &data) < 0)
5671 goto error;
5672 if (pick_side(local, &data) < 0)
5673 goto error;
5675 local->side++;
5676 level++;
5677 init = 1;
5680 clear_lexmin_data(&data);
5681 isl_basic_set_free(bset);
5683 return data.sol;
5684 error:
5685 clear_lexmin_data(&data);
5686 isl_basic_set_free(bset);
5687 isl_vec_free(data.sol);
5688 return NULL;
5691 /* Wrapper for a tableau that is used for computing
5692 * the lexicographically smallest rational point of a non-negative set.
5693 * This point is represented by the sample value of "tab",
5694 * unless "tab" is empty.
5696 struct isl_tab_lexmin {
5697 isl_ctx *ctx;
5698 struct isl_tab *tab;
5701 /* Free "tl" and return NULL.
5703 __isl_null isl_tab_lexmin *isl_tab_lexmin_free(__isl_take isl_tab_lexmin *tl)
5705 if (!tl)
5706 return NULL;
5707 isl_ctx_deref(tl->ctx);
5708 isl_tab_free(tl->tab);
5709 free(tl);
5711 return NULL;
5714 /* Construct an isl_tab_lexmin for computing
5715 * the lexicographically smallest rational point in "bset",
5716 * assuming that all variables are non-negative.
5718 __isl_give isl_tab_lexmin *isl_tab_lexmin_from_basic_set(
5719 __isl_take isl_basic_set *bset)
5721 isl_ctx *ctx;
5722 isl_tab_lexmin *tl;
5724 if (!bset)
5725 return NULL;
5727 ctx = isl_basic_set_get_ctx(bset);
5728 tl = isl_calloc_type(ctx, struct isl_tab_lexmin);
5729 if (!tl)
5730 goto error;
5731 tl->ctx = ctx;
5732 isl_ctx_ref(ctx);
5733 tl->tab = tab_for_lexmin(bset, NULL, 0, 0);
5734 isl_basic_set_free(bset);
5735 if (!tl->tab)
5736 return isl_tab_lexmin_free(tl);
5737 return tl;
5738 error:
5739 isl_basic_set_free(bset);
5740 isl_tab_lexmin_free(tl);
5741 return NULL;
5744 /* Return the dimension of the set represented by "tl".
5746 int isl_tab_lexmin_dim(__isl_keep isl_tab_lexmin *tl)
5748 return tl ? tl->tab->n_var : -1;
5751 /* Add the equality with coefficients "eq" to "tl", updating the optimal
5752 * solution if needed.
5753 * The equality is added as two opposite inequality constraints.
5755 __isl_give isl_tab_lexmin *isl_tab_lexmin_add_eq(__isl_take isl_tab_lexmin *tl,
5756 isl_int *eq)
5758 unsigned n_var;
5760 if (!tl || !eq)
5761 return isl_tab_lexmin_free(tl);
5763 if (isl_tab_extend_cons(tl->tab, 2) < 0)
5764 return isl_tab_lexmin_free(tl);
5765 n_var = tl->tab->n_var;
5766 isl_seq_neg(eq, eq, 1 + n_var);
5767 tl->tab = add_lexmin_ineq(tl->tab, eq);
5768 isl_seq_neg(eq, eq, 1 + n_var);
5769 tl->tab = add_lexmin_ineq(tl->tab, eq);
5771 if (!tl->tab)
5772 return isl_tab_lexmin_free(tl);
5774 return tl;
5777 /* Add cuts to "tl" until the sample value reaches an integer value or
5778 * until the result becomes empty.
5780 __isl_give isl_tab_lexmin *isl_tab_lexmin_cut_to_integer(
5781 __isl_take isl_tab_lexmin *tl)
5783 if (!tl)
5784 return NULL;
5785 tl->tab = cut_to_integer_lexmin(tl->tab, CUT_ONE);
5786 if (!tl->tab)
5787 return isl_tab_lexmin_free(tl);
5788 return tl;
5791 /* Return the lexicographically smallest rational point in the basic set
5792 * from which "tl" was constructed.
5793 * If the original input was empty, then return a zero-length vector.
5795 __isl_give isl_vec *isl_tab_lexmin_get_solution(__isl_keep isl_tab_lexmin *tl)
5797 if (!tl)
5798 return NULL;
5799 if (tl->tab->empty)
5800 return isl_vec_alloc(tl->ctx, 0);
5801 else
5802 return isl_tab_get_sample_value(tl->tab);
5805 struct isl_sol_pma {
5806 struct isl_sol sol;
5807 isl_pw_multi_aff *pma;
5808 isl_set *empty;
5811 static void sol_pma_free(struct isl_sol *sol)
5813 struct isl_sol_pma *sol_pma = (struct isl_sol_pma *) sol;
5814 isl_pw_multi_aff_free(sol_pma->pma);
5815 isl_set_free(sol_pma->empty);
5818 /* This function is called for parts of the context where there is
5819 * no solution, with "bset" corresponding to the context tableau.
5820 * Simply add the basic set to the set "empty".
5822 static void sol_pma_add_empty(struct isl_sol_pma *sol,
5823 __isl_take isl_basic_set *bset)
5825 if (!bset || !sol->empty)
5826 goto error;
5828 sol->empty = isl_set_grow(sol->empty, 1);
5829 bset = isl_basic_set_simplify(bset);
5830 bset = isl_basic_set_finalize(bset);
5831 sol->empty = isl_set_add_basic_set(sol->empty, bset);
5832 if (!sol->empty)
5833 sol->sol.error = 1;
5834 return;
5835 error:
5836 isl_basic_set_free(bset);
5837 sol->sol.error = 1;
5840 /* Given a basic set "dom" that represents the context and a tuple of
5841 * affine expressions "maff" defined over this domain, construct
5842 * an isl_pw_multi_aff with a single cell corresponding to "dom" and
5843 * the affine expressions in "maff".
5845 static void sol_pma_add(struct isl_sol_pma *sol,
5846 __isl_take isl_basic_set *dom, __isl_take isl_multi_aff *maff)
5848 isl_pw_multi_aff *pma;
5850 dom = isl_basic_set_simplify(dom);
5851 dom = isl_basic_set_finalize(dom);
5852 pma = isl_pw_multi_aff_alloc(isl_set_from_basic_set(dom), maff);
5853 sol->pma = isl_pw_multi_aff_add_disjoint(sol->pma, pma);
5854 if (!sol->pma)
5855 sol->sol.error = 1;
5858 static void sol_pma_add_empty_wrap(struct isl_sol *sol,
5859 __isl_take isl_basic_set *bset)
5861 sol_pma_add_empty((struct isl_sol_pma *)sol, bset);
5864 static void sol_pma_add_wrap(struct isl_sol *sol,
5865 __isl_take isl_basic_set *dom, __isl_take isl_multi_aff *ma)
5867 sol_pma_add((struct isl_sol_pma *)sol, dom, ma);
5870 /* Construct an isl_sol_pma structure for accumulating the solution.
5871 * If track_empty is set, then we also keep track of the parts
5872 * of the context where there is no solution.
5873 * If max is set, then we are solving a maximization, rather than
5874 * a minimization problem, which means that the variables in the
5875 * tableau have value "M - x" rather than "M + x".
5877 static struct isl_sol *sol_pma_init(__isl_keep isl_basic_map *bmap,
5878 __isl_take isl_basic_set *dom, int track_empty, int max)
5880 struct isl_sol_pma *sol_pma = NULL;
5881 isl_space *space;
5883 if (!bmap)
5884 goto error;
5886 sol_pma = isl_calloc_type(bmap->ctx, struct isl_sol_pma);
5887 if (!sol_pma)
5888 goto error;
5890 sol_pma->sol.free = &sol_pma_free;
5891 if (sol_init(&sol_pma->sol, bmap, dom, max) < 0)
5892 goto error;
5893 sol_pma->sol.add = &sol_pma_add_wrap;
5894 sol_pma->sol.add_empty = track_empty ? &sol_pma_add_empty_wrap : NULL;
5895 space = isl_space_copy(sol_pma->sol.space);
5896 sol_pma->pma = isl_pw_multi_aff_empty(space);
5897 if (!sol_pma->pma)
5898 goto error;
5900 if (track_empty) {
5901 sol_pma->empty = isl_set_alloc_space(isl_basic_set_get_space(dom),
5902 1, ISL_SET_DISJOINT);
5903 if (!sol_pma->empty)
5904 goto error;
5907 isl_basic_set_free(dom);
5908 return &sol_pma->sol;
5909 error:
5910 isl_basic_set_free(dom);
5911 sol_free(&sol_pma->sol);
5912 return NULL;
5915 /* Base case of isl_tab_basic_map_partial_lexopt, after removing
5916 * some obvious symmetries.
5918 * We call basic_map_partial_lexopt_base_sol and extract the results.
5920 static __isl_give isl_pw_multi_aff *basic_map_partial_lexopt_base_pw_multi_aff(
5921 __isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom,
5922 __isl_give isl_set **empty, int max)
5924 isl_pw_multi_aff *result = NULL;
5925 struct isl_sol *sol;
5926 struct isl_sol_pma *sol_pma;
5928 sol = basic_map_partial_lexopt_base_sol(bmap, dom, empty, max,
5929 &sol_pma_init);
5930 if (!sol)
5931 return NULL;
5932 sol_pma = (struct isl_sol_pma *) sol;
5934 result = isl_pw_multi_aff_copy(sol_pma->pma);
5935 if (empty)
5936 *empty = isl_set_copy(sol_pma->empty);
5937 sol_free(&sol_pma->sol);
5938 return result;
5941 /* Given that the last input variable of "maff" represents the minimum
5942 * of some bounds, check whether we need to plug in the expression
5943 * of the minimum.
5945 * In particular, check if the last input variable appears in any
5946 * of the expressions in "maff".
5948 static isl_bool need_substitution(__isl_keep isl_multi_aff *maff)
5950 int i;
5951 isl_size n_in;
5952 unsigned pos;
5954 n_in = isl_multi_aff_dim(maff, isl_dim_in);
5955 if (n_in < 0)
5956 return isl_bool_error;
5957 pos = n_in - 1;
5959 for (i = 0; i < maff->n; ++i) {
5960 isl_bool involves;
5962 involves = isl_aff_involves_dims(maff->u.p[i],
5963 isl_dim_in, pos, 1);
5964 if (involves < 0 || involves)
5965 return involves;
5968 return isl_bool_false;
5971 /* Given a set of upper bounds on the last "input" variable m,
5972 * construct a piecewise affine expression that selects
5973 * the minimal upper bound to m, i.e.,
5974 * divide the space into cells where one
5975 * of the upper bounds is smaller than all the others and select
5976 * this upper bound on that cell.
5978 * In particular, if there are n bounds b_i, then the result
5979 * consists of n cell, each one of the form
5981 * b_i <= b_j for j > i
5982 * b_i < b_j for j < i
5984 * The affine expression on this cell is
5986 * b_i
5988 static __isl_give isl_pw_aff *set_minimum_pa(__isl_take isl_space *space,
5989 __isl_take isl_mat *var)
5991 int i;
5992 isl_aff *aff = NULL;
5993 isl_basic_set *bset = NULL;
5994 isl_pw_aff *paff = NULL;
5995 isl_space *pw_space;
5996 isl_local_space *ls = NULL;
5998 if (!space || !var)
5999 goto error;
6001 ls = isl_local_space_from_space(isl_space_copy(space));
6002 pw_space = isl_space_copy(space);
6003 pw_space = isl_space_from_domain(pw_space);
6004 pw_space = isl_space_add_dims(pw_space, isl_dim_out, 1);
6005 paff = isl_pw_aff_alloc_size(pw_space, var->n_row);
6007 for (i = 0; i < var->n_row; ++i) {
6008 isl_pw_aff *paff_i;
6010 aff = isl_aff_alloc(isl_local_space_copy(ls));
6011 bset = isl_basic_set_alloc_space(isl_space_copy(space), 0,
6012 0, var->n_row - 1);
6013 if (!aff || !bset)
6014 goto error;
6015 isl_int_set_si(aff->v->el[0], 1);
6016 isl_seq_cpy(aff->v->el + 1, var->row[i], var->n_col);
6017 isl_int_set_si(aff->v->el[1 + var->n_col], 0);
6018 bset = select_minimum(bset, var, i);
6019 paff_i = isl_pw_aff_alloc(isl_set_from_basic_set(bset), aff);
6020 paff = isl_pw_aff_add_disjoint(paff, paff_i);
6023 isl_local_space_free(ls);
6024 isl_space_free(space);
6025 isl_mat_free(var);
6026 return paff;
6027 error:
6028 isl_aff_free(aff);
6029 isl_basic_set_free(bset);
6030 isl_pw_aff_free(paff);
6031 isl_local_space_free(ls);
6032 isl_space_free(space);
6033 isl_mat_free(var);
6034 return NULL;
6037 /* Given a piecewise multi-affine expression of which the last input variable
6038 * is the minimum of the bounds in "cst", plug in the value of the minimum.
6039 * This minimum expression is given in "min_expr_pa".
6040 * The set "min_expr" contains the same information, but in the form of a set.
6041 * The variable is subsequently projected out.
6043 * The implementation is similar to those of "split" and "split_domain".
6044 * If the variable appears in a given expression, then minimum expression
6045 * is plugged in. Otherwise, if the variable appears in the constraints
6046 * and a split is required, then the domain is split. Otherwise, no split
6047 * is performed.
6049 static __isl_give isl_pw_multi_aff *split_domain_pma(
6050 __isl_take isl_pw_multi_aff *opt, __isl_take isl_pw_aff *min_expr_pa,
6051 __isl_take isl_set *min_expr, __isl_take isl_mat *cst)
6053 isl_size n_in;
6054 int i;
6055 isl_space *space;
6056 isl_pw_multi_aff *res;
6058 if (!opt || !min_expr || !cst)
6059 goto error;
6061 n_in = isl_pw_multi_aff_dim(opt, isl_dim_in);
6062 if (n_in < 0)
6063 goto error;
6064 space = isl_pw_multi_aff_get_space(opt);
6065 space = isl_space_drop_dims(space, isl_dim_in, n_in - 1, 1);
6066 res = isl_pw_multi_aff_empty(space);
6068 for (i = 0; i < opt->n; ++i) {
6069 isl_bool subs;
6070 isl_pw_multi_aff *pma;
6072 pma = isl_pw_multi_aff_alloc(isl_set_copy(opt->p[i].set),
6073 isl_multi_aff_copy(opt->p[i].maff));
6074 subs = need_substitution(opt->p[i].maff);
6075 if (subs < 0) {
6076 pma = isl_pw_multi_aff_free(pma);
6077 } else if (subs) {
6078 pma = isl_pw_multi_aff_substitute(pma,
6079 n_in - 1, min_expr_pa);
6080 } else {
6081 isl_bool split;
6082 split = need_split_set(opt->p[i].set, cst);
6083 if (split < 0)
6084 pma = isl_pw_multi_aff_free(pma);
6085 else if (split)
6086 pma = isl_pw_multi_aff_intersect_domain(pma,
6087 isl_set_copy(min_expr));
6089 pma = isl_pw_multi_aff_project_out(pma,
6090 isl_dim_in, n_in - 1, 1);
6092 res = isl_pw_multi_aff_add_disjoint(res, pma);
6095 isl_pw_multi_aff_free(opt);
6096 isl_pw_aff_free(min_expr_pa);
6097 isl_set_free(min_expr);
6098 isl_mat_free(cst);
6099 return res;
6100 error:
6101 isl_pw_multi_aff_free(opt);
6102 isl_pw_aff_free(min_expr_pa);
6103 isl_set_free(min_expr);
6104 isl_mat_free(cst);
6105 return NULL;
6108 static __isl_give isl_pw_multi_aff *basic_map_partial_lexopt_pw_multi_aff(
6109 __isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom,
6110 __isl_give isl_set **empty, int max);
6112 /* This function is called from basic_map_partial_lexopt_symm.
6113 * The last variable of "bmap" and "dom" corresponds to the minimum
6114 * of the bounds in "cst". "map_space" is the space of the original
6115 * input relation (of basic_map_partial_lexopt_symm) and "set_space"
6116 * is the space of the original domain.
6118 * We recursively call basic_map_partial_lexopt and then plug in
6119 * the definition of the minimum in the result.
6121 static __isl_give isl_pw_multi_aff *
6122 basic_map_partial_lexopt_symm_core_pw_multi_aff(
6123 __isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom,
6124 __isl_give isl_set **empty, int max, __isl_take isl_mat *cst,
6125 __isl_take isl_space *map_space, __isl_take isl_space *set_space)
6127 isl_pw_multi_aff *opt;
6128 isl_pw_aff *min_expr_pa;
6129 isl_set *min_expr;
6131 min_expr = set_minimum(isl_basic_set_get_space(dom), isl_mat_copy(cst));
6132 min_expr_pa = set_minimum_pa(isl_basic_set_get_space(dom),
6133 isl_mat_copy(cst));
6135 opt = basic_map_partial_lexopt_pw_multi_aff(bmap, dom, empty, max);
6137 if (empty) {
6138 *empty = split(*empty,
6139 isl_set_copy(min_expr), isl_mat_copy(cst));
6140 *empty = isl_set_reset_space(*empty, set_space);
6143 opt = split_domain_pma(opt, min_expr_pa, min_expr, cst);
6144 opt = isl_pw_multi_aff_reset_space(opt, map_space);
6146 return opt;
6149 #undef TYPE
6150 #define TYPE isl_pw_multi_aff
6151 #undef SUFFIX
6152 #define SUFFIX _pw_multi_aff
6153 #include "isl_tab_lexopt_templ.c"