scheduler.c: use isl_bool enum instead of plain integers
[isl.git] / isl_tab_pip.c
blobe305cfd933274ab46829919c70ec1c0793992889
1 /*
2 * Copyright 2008-2009 Katholieke Universiteit Leuven
3 * Copyright 2010 INRIA Saclay
4 * Copyright 2016-2017 Sven Verdoolaege
6 * Use of this software is governed by the MIT license
8 * Written by Sven Verdoolaege, K.U.Leuven, Departement
9 * Computerwetenschappen, Celestijnenlaan 200A, B-3001 Leuven, Belgium
10 * and INRIA Saclay - Ile-de-France, Parc Club Orsay Universite,
11 * ZAC des vignes, 4 rue Jacques Monod, 91893 Orsay, France
14 #include <isl_ctx_private.h>
15 #include "isl_map_private.h"
16 #include <isl_seq.h>
17 #include "isl_tab.h"
18 #include "isl_sample.h"
19 #include <isl_mat_private.h>
20 #include <isl_vec_private.h>
21 #include <isl_aff_private.h>
22 #include <isl_constraint_private.h>
23 #include <isl_options_private.h>
24 #include <isl_config.h>
26 #include <bset_to_bmap.c>
29 * The implementation of parametric integer linear programming in this file
30 * was inspired by the paper "Parametric Integer Programming" and the
31 * report "Solving systems of affine (in)equalities" by Paul Feautrier
32 * (and others).
34 * The strategy used for obtaining a feasible solution is different
35 * from the one used in isl_tab.c. In particular, in isl_tab.c,
36 * upon finding a constraint that is not yet satisfied, we pivot
37 * in a row that increases the constant term of the row holding the
38 * constraint, making sure the sample solution remains feasible
39 * for all the constraints it already satisfied.
40 * Here, we always pivot in the row holding the constraint,
41 * choosing a column that induces the lexicographically smallest
42 * increment to the sample solution.
44 * By starting out from a sample value that is lexicographically
45 * smaller than any integer point in the problem space, the first
46 * feasible integer sample point we find will also be the lexicographically
47 * smallest. If all variables can be assumed to be non-negative,
48 * then the initial sample value may be chosen equal to zero.
49 * However, we will not make this assumption. Instead, we apply
50 * the "big parameter" trick. Any variable x is then not directly
51 * used in the tableau, but instead it is represented by another
52 * variable x' = M + x, where M is an arbitrarily large (positive)
53 * value. x' is therefore always non-negative, whatever the value of x.
54 * Taking as initial sample value x' = 0 corresponds to x = -M,
55 * which is always smaller than any possible value of x.
57 * The big parameter trick is used in the main tableau and
58 * also in the context tableau if isl_context_lex is used.
59 * In this case, each tableaus has its own big parameter.
60 * Before doing any real work, we check if all the parameters
61 * happen to be non-negative. If so, we drop the column corresponding
62 * to M from the initial context tableau.
63 * If isl_context_gbr is used, then the big parameter trick is only
64 * used in the main tableau.
67 struct isl_context;
68 struct isl_context_op {
69 /* detect nonnegative parameters in context and mark them in tab */
70 struct isl_tab *(*detect_nonnegative_parameters)(
71 struct isl_context *context, struct isl_tab *tab);
72 /* return temporary reference to basic set representation of context */
73 struct isl_basic_set *(*peek_basic_set)(struct isl_context *context);
74 /* return temporary reference to tableau representation of context */
75 struct isl_tab *(*peek_tab)(struct isl_context *context);
76 /* add equality; check is 1 if eq may not be valid;
77 * update is 1 if we may want to call ineq_sign on context later.
79 void (*add_eq)(struct isl_context *context, isl_int *eq,
80 int check, int update);
81 /* add inequality; check is 1 if ineq may not be valid;
82 * update is 1 if we may want to call ineq_sign on context later.
84 void (*add_ineq)(struct isl_context *context, isl_int *ineq,
85 int check, int update);
86 /* check sign of ineq based on previous information.
87 * strict is 1 if saturation should be treated as a positive sign.
89 enum isl_tab_row_sign (*ineq_sign)(struct isl_context *context,
90 isl_int *ineq, int strict);
91 /* check if inequality maintains feasibility */
92 int (*test_ineq)(struct isl_context *context, isl_int *ineq);
93 /* return index of a div that corresponds to "div" */
94 int (*get_div)(struct isl_context *context, struct isl_tab *tab,
95 struct isl_vec *div);
96 /* insert div "div" to context at "pos" and return non-negativity */
97 isl_bool (*insert_div)(struct isl_context *context, int pos,
98 __isl_keep isl_vec *div);
99 int (*detect_equalities)(struct isl_context *context,
100 struct isl_tab *tab);
101 /* return row index of "best" split */
102 int (*best_split)(struct isl_context *context, struct isl_tab *tab);
103 /* check if context has already been determined to be empty */
104 int (*is_empty)(struct isl_context *context);
105 /* check if context is still usable */
106 int (*is_ok)(struct isl_context *context);
107 /* save a copy/snapshot of context */
108 void *(*save)(struct isl_context *context);
109 /* restore saved context */
110 void (*restore)(struct isl_context *context, void *);
111 /* discard saved context */
112 void (*discard)(void *);
113 /* invalidate context */
114 void (*invalidate)(struct isl_context *context);
115 /* free context */
116 __isl_null struct isl_context *(*free)(struct isl_context *context);
119 /* Shared parts of context representation.
121 * "n_unknown" is the number of final unknown integer divisions
122 * in the input domain.
124 struct isl_context {
125 struct isl_context_op *op;
126 int n_unknown;
129 struct isl_context_lex {
130 struct isl_context context;
131 struct isl_tab *tab;
134 /* A stack (linked list) of solutions of subtrees of the search space.
136 * "ma" describes the solution as a function of "dom".
137 * In particular, the domain space of "ma" is equal to the space of "dom".
139 * If "ma" is NULL, then there is no solution on "dom".
141 struct isl_partial_sol {
142 int level;
143 struct isl_basic_set *dom;
144 isl_multi_aff *ma;
146 struct isl_partial_sol *next;
149 struct isl_sol;
150 struct isl_sol_callback {
151 struct isl_tab_callback callback;
152 struct isl_sol *sol;
155 /* isl_sol is an interface for constructing a solution to
156 * a parametric integer linear programming problem.
157 * Every time the algorithm reaches a state where a solution
158 * can be read off from the tableau, the function "add" is called
159 * on the isl_sol passed to find_solutions_main. In a state where
160 * the tableau is empty, "add_empty" is called instead.
161 * "free" is called to free the implementation specific fields, if any.
163 * "error" is set if some error has occurred. This flag invalidates
164 * the remainder of the data structure.
165 * If "rational" is set, then a rational optimization is being performed.
166 * "level" is the current level in the tree with nodes for each
167 * split in the context.
168 * If "max" is set, then a maximization problem is being solved, rather than
169 * a minimization problem, which means that the variables in the
170 * tableau have value "M - x" rather than "M + x".
171 * "n_out" is the number of output dimensions in the input.
172 * "space" is the space in which the solution (and also the input) lives.
174 * The context tableau is owned by isl_sol and is updated incrementally.
176 * There are currently two implementations of this interface,
177 * isl_sol_map, which simply collects the solutions in an isl_map
178 * and (optionally) the parts of the context where there is no solution
179 * in an isl_set, and
180 * isl_sol_pma, which collects an isl_pw_multi_aff instead.
182 struct isl_sol {
183 int error;
184 int rational;
185 int level;
186 int max;
187 int n_out;
188 isl_space *space;
189 struct isl_context *context;
190 struct isl_partial_sol *partial;
191 void (*add)(struct isl_sol *sol,
192 __isl_take isl_basic_set *dom, __isl_take isl_multi_aff *ma);
193 void (*add_empty)(struct isl_sol *sol, struct isl_basic_set *bset);
194 void (*free)(struct isl_sol *sol);
195 struct isl_sol_callback dec_level;
198 static void sol_free(struct isl_sol *sol)
200 struct isl_partial_sol *partial, *next;
201 if (!sol)
202 return;
203 for (partial = sol->partial; partial; partial = next) {
204 next = partial->next;
205 isl_basic_set_free(partial->dom);
206 isl_multi_aff_free(partial->ma);
207 free(partial);
209 isl_space_free(sol->space);
210 if (sol->context)
211 sol->context->op->free(sol->context);
212 sol->free(sol);
213 free(sol);
216 /* Push a partial solution represented by a domain and function "ma"
217 * onto the stack of partial solutions.
218 * If "ma" is NULL, then "dom" represents a part of the domain
219 * with no solution.
221 static void sol_push_sol(struct isl_sol *sol,
222 __isl_take isl_basic_set *dom, __isl_take isl_multi_aff *ma)
224 struct isl_partial_sol *partial;
226 if (sol->error || !dom)
227 goto error;
229 partial = isl_alloc_type(dom->ctx, struct isl_partial_sol);
230 if (!partial)
231 goto error;
233 partial->level = sol->level;
234 partial->dom = dom;
235 partial->ma = ma;
236 partial->next = sol->partial;
238 sol->partial = partial;
240 return;
241 error:
242 isl_basic_set_free(dom);
243 isl_multi_aff_free(ma);
244 sol->error = 1;
247 /* Check that the final columns of "M", starting at "first", are zero.
249 static isl_stat check_final_columns_are_zero(__isl_keep isl_mat *M,
250 unsigned first)
252 int i;
253 unsigned rows, cols, n;
255 if (!M)
256 return isl_stat_error;
257 rows = isl_mat_rows(M);
258 cols = isl_mat_cols(M);
259 n = cols - first;
260 for (i = 0; i < rows; ++i)
261 if (isl_seq_first_non_zero(M->row[i] + first, n) != -1)
262 isl_die(isl_mat_get_ctx(M), isl_error_internal,
263 "final columns should be zero",
264 return isl_stat_error);
265 return isl_stat_ok;
268 /* Set the affine expressions in "ma" according to the rows in "M", which
269 * are defined over the local space "ls".
270 * The matrix "M" may have extra (zero) columns beyond the number
271 * of variables in "ls".
273 static __isl_give isl_multi_aff *set_from_affine_matrix(
274 __isl_take isl_multi_aff *ma, __isl_take isl_local_space *ls,
275 __isl_take isl_mat *M)
277 int i, dim;
278 isl_aff *aff;
280 if (!ma || !ls || !M)
281 goto error;
283 dim = isl_local_space_dim(ls, isl_dim_all);
284 if (check_final_columns_are_zero(M, 1 + dim) < 0)
285 goto error;
286 for (i = 1; i < M->n_row; ++i) {
287 aff = isl_aff_alloc(isl_local_space_copy(ls));
288 if (aff) {
289 isl_int_set(aff->v->el[0], M->row[0][0]);
290 isl_seq_cpy(aff->v->el + 1, M->row[i], 1 + dim);
292 aff = isl_aff_normalize(aff);
293 ma = isl_multi_aff_set_aff(ma, i - 1, aff);
295 isl_local_space_free(ls);
296 isl_mat_free(M);
298 return ma;
299 error:
300 isl_local_space_free(ls);
301 isl_mat_free(M);
302 isl_multi_aff_free(ma);
303 return NULL;
306 /* Push a partial solution represented by a domain and mapping M
307 * onto the stack of partial solutions.
309 * The affine matrix "M" maps the dimensions of the context
310 * to the output variables. Convert it into an isl_multi_aff and
311 * then call sol_push_sol.
313 * Note that the description of the initial context may have involved
314 * existentially quantified variables, in which case they also appear
315 * in "dom". These need to be removed before creating the affine
316 * expression because an affine expression cannot be defined in terms
317 * of existentially quantified variables without a known representation.
318 * Since newly added integer divisions are inserted before these
319 * existentially quantified variables, they are still in the final
320 * positions and the corresponding final columns of "M" are zero
321 * because align_context_divs adds the existentially quantified
322 * variables of the context to the main tableau without any constraints and
323 * any equality constraints that are added later on can only serve
324 * to eliminate these existentially quantified variables.
326 static void sol_push_sol_mat(struct isl_sol *sol,
327 __isl_take isl_basic_set *dom, __isl_take isl_mat *M)
329 isl_local_space *ls;
330 isl_multi_aff *ma;
331 int n_div, n_known;
333 n_div = isl_basic_set_dim(dom, isl_dim_div);
334 n_known = n_div - sol->context->n_unknown;
336 ma = isl_multi_aff_alloc(isl_space_copy(sol->space));
337 ls = isl_basic_set_get_local_space(dom);
338 ls = isl_local_space_drop_dims(ls, isl_dim_div,
339 n_known, n_div - n_known);
340 ma = set_from_affine_matrix(ma, ls, M);
342 if (!ma)
343 dom = isl_basic_set_free(dom);
344 sol_push_sol(sol, dom, ma);
347 /* Pop one partial solution from the partial solution stack and
348 * pass it on to sol->add or sol->add_empty.
350 static void sol_pop_one(struct isl_sol *sol)
352 struct isl_partial_sol *partial;
354 partial = sol->partial;
355 sol->partial = partial->next;
357 if (partial->ma)
358 sol->add(sol, partial->dom, partial->ma);
359 else
360 sol->add_empty(sol, partial->dom);
361 free(partial);
364 /* Return a fresh copy of the domain represented by the context tableau.
366 static struct isl_basic_set *sol_domain(struct isl_sol *sol)
368 struct isl_basic_set *bset;
370 if (sol->error)
371 return NULL;
373 bset = isl_basic_set_dup(sol->context->op->peek_basic_set(sol->context));
374 bset = isl_basic_set_update_from_tab(bset,
375 sol->context->op->peek_tab(sol->context));
377 return bset;
380 /* Check whether two partial solutions have the same affine expressions.
382 static isl_bool same_solution(struct isl_partial_sol *s1,
383 struct isl_partial_sol *s2)
385 if (!s1->ma != !s2->ma)
386 return isl_bool_false;
387 if (!s1->ma)
388 return isl_bool_true;
390 return isl_multi_aff_plain_is_equal(s1->ma, s2->ma);
393 /* Swap the initial two partial solutions in "sol".
395 * That is, go from
397 * sol->partial = p1; p1->next = p2; p2->next = p3
399 * to
401 * sol->partial = p2; p2->next = p1; p1->next = p3
403 static void swap_initial(struct isl_sol *sol)
405 struct isl_partial_sol *partial;
407 partial = sol->partial;
408 sol->partial = partial->next;
409 partial->next = partial->next->next;
410 sol->partial->next = partial;
413 /* Combine the initial two partial solution of "sol" into
414 * a partial solution with the current context domain of "sol" and
415 * the function description of the second partial solution in the list.
416 * The level of the new partial solution is set to the current level.
418 * That is, the first two partial solutions (D1,M1) and (D2,M2) are
419 * replaced by (D,M2), where D is the domain of "sol", which is assumed
420 * to be the union of D1 and D2, while M1 is assumed to be equal to M2
421 * (at least on D1).
423 static isl_stat combine_initial_into_second(struct isl_sol *sol)
425 struct isl_partial_sol *partial;
426 isl_basic_set *bset;
428 partial = sol->partial;
430 bset = sol_domain(sol);
431 isl_basic_set_free(partial->next->dom);
432 partial->next->dom = bset;
433 partial->next->level = sol->level;
435 if (!bset)
436 return isl_stat_error;
438 sol->partial = partial->next;
439 isl_basic_set_free(partial->dom);
440 isl_multi_aff_free(partial->ma);
441 free(partial);
443 return isl_stat_ok;
446 /* Are "ma1" and "ma2" equal to each other on "dom"?
448 * Combine "ma1" and "ma2" with "dom" and check if the results are the same.
449 * "dom" may have existentially quantified variables. Eliminate them first
450 * as otherwise they would have to be eliminated twice, in a more complicated
451 * context.
453 static isl_bool equal_on_domain(__isl_keep isl_multi_aff *ma1,
454 __isl_keep isl_multi_aff *ma2, __isl_keep isl_basic_set *dom)
456 isl_set *set;
457 isl_pw_multi_aff *pma1, *pma2;
458 isl_bool equal;
460 set = isl_basic_set_compute_divs(isl_basic_set_copy(dom));
461 pma1 = isl_pw_multi_aff_alloc(isl_set_copy(set),
462 isl_multi_aff_copy(ma1));
463 pma2 = isl_pw_multi_aff_alloc(set, isl_multi_aff_copy(ma2));
464 equal = isl_pw_multi_aff_is_equal(pma1, pma2);
465 isl_pw_multi_aff_free(pma1);
466 isl_pw_multi_aff_free(pma2);
468 return equal;
471 /* The initial two partial solutions of "sol" are known to be at
472 * the same level.
473 * If they represent the same solution (on different parts of the domain),
474 * then combine them into a single solution at the current level.
475 * Otherwise, pop them both.
477 * Even if the two partial solution are not obviously the same,
478 * one may still be a simplification of the other over its own domain.
479 * Also check if the two sets of affine functions are equal when
480 * restricted to one of the domains. If so, combine the two
481 * using the set of affine functions on the other domain.
482 * That is, for two partial solutions (D1,M1) and (D2,M2),
483 * if M1 = M2 on D1, then the pair of partial solutions can
484 * be replaced by (D1+D2,M2) and similarly when M1 = M2 on D2.
486 static isl_stat combine_initial_if_equal(struct isl_sol *sol)
488 struct isl_partial_sol *partial;
489 isl_bool same;
491 partial = sol->partial;
493 same = same_solution(partial, partial->next);
494 if (same < 0)
495 return isl_stat_error;
496 if (same)
497 return combine_initial_into_second(sol);
498 if (partial->ma && partial->next->ma) {
499 same = equal_on_domain(partial->ma, partial->next->ma,
500 partial->dom);
501 if (same < 0)
502 return isl_stat_error;
503 if (same)
504 return combine_initial_into_second(sol);
505 same = equal_on_domain(partial->ma, partial->next->ma,
506 partial->next->dom);
507 if (same) {
508 swap_initial(sol);
509 return combine_initial_into_second(sol);
513 sol_pop_one(sol);
514 sol_pop_one(sol);
516 return isl_stat_ok;
519 /* Pop all solutions from the partial solution stack that were pushed onto
520 * the stack at levels that are deeper than the current level.
521 * If the two topmost elements on the stack have the same level
522 * and represent the same solution, then their domains are combined.
523 * This combined domain is the same as the current context domain
524 * as sol_pop is called each time we move back to a higher level.
525 * If the outer level (0) has been reached, then all partial solutions
526 * at the current level are also popped off.
528 static void sol_pop(struct isl_sol *sol)
530 struct isl_partial_sol *partial;
532 if (sol->error)
533 return;
535 partial = sol->partial;
536 if (!partial)
537 return;
539 if (partial->level == 0 && sol->level == 0) {
540 for (partial = sol->partial; partial; partial = sol->partial)
541 sol_pop_one(sol);
542 return;
545 if (partial->level <= sol->level)
546 return;
548 if (partial->next && partial->next->level == partial->level) {
549 if (combine_initial_if_equal(sol) < 0)
550 goto error;
551 } else
552 sol_pop_one(sol);
554 if (sol->level == 0) {
555 for (partial = sol->partial; partial; partial = sol->partial)
556 sol_pop_one(sol);
557 return;
560 if (0)
561 error: sol->error = 1;
564 static void sol_dec_level(struct isl_sol *sol)
566 if (sol->error)
567 return;
569 sol->level--;
571 sol_pop(sol);
574 static isl_stat sol_dec_level_wrap(struct isl_tab_callback *cb)
576 struct isl_sol_callback *callback = (struct isl_sol_callback *)cb;
578 sol_dec_level(callback->sol);
580 return callback->sol->error ? isl_stat_error : isl_stat_ok;
583 /* Move down to next level and push callback onto context tableau
584 * to decrease the level again when it gets rolled back across
585 * the current state. That is, dec_level will be called with
586 * the context tableau in the same state as it is when inc_level
587 * is called.
589 static void sol_inc_level(struct isl_sol *sol)
591 struct isl_tab *tab;
593 if (sol->error)
594 return;
596 sol->level++;
597 tab = sol->context->op->peek_tab(sol->context);
598 if (isl_tab_push_callback(tab, &sol->dec_level.callback) < 0)
599 sol->error = 1;
602 static void scale_rows(struct isl_mat *mat, isl_int m, int n_row)
604 int i;
606 if (isl_int_is_one(m))
607 return;
609 for (i = 0; i < n_row; ++i)
610 isl_seq_scale(mat->row[i], mat->row[i], m, mat->n_col);
613 /* Add the solution identified by the tableau and the context tableau.
615 * The layout of the variables is as follows.
616 * tab->n_var is equal to the total number of variables in the input
617 * map (including divs that were copied from the context)
618 * + the number of extra divs constructed
619 * Of these, the first tab->n_param and the last tab->n_div variables
620 * correspond to the variables in the context, i.e.,
621 * tab->n_param + tab->n_div = context_tab->n_var
622 * tab->n_param is equal to the number of parameters and input
623 * dimensions in the input map
624 * tab->n_div is equal to the number of divs in the context
626 * If there is no solution, then call add_empty with a basic set
627 * that corresponds to the context tableau. (If add_empty is NULL,
628 * then do nothing).
630 * If there is a solution, then first construct a matrix that maps
631 * all dimensions of the context to the output variables, i.e.,
632 * the output dimensions in the input map.
633 * The divs in the input map (if any) that do not correspond to any
634 * div in the context do not appear in the solution.
635 * The algorithm will make sure that they have an integer value,
636 * but these values themselves are of no interest.
637 * We have to be careful not to drop or rearrange any divs in the
638 * context because that would change the meaning of the matrix.
640 * To extract the value of the output variables, it should be noted
641 * that we always use a big parameter M in the main tableau and so
642 * the variable stored in this tableau is not an output variable x itself, but
643 * x' = M + x (in case of minimization)
644 * or
645 * x' = M - x (in case of maximization)
646 * If x' appears in a column, then its optimal value is zero,
647 * which means that the optimal value of x is an unbounded number
648 * (-M for minimization and M for maximization).
649 * We currently assume that the output dimensions in the original map
650 * are bounded, so this cannot occur.
651 * Similarly, when x' appears in a row, then the coefficient of M in that
652 * row is necessarily 1.
653 * If the row in the tableau represents
654 * d x' = c + d M + e(y)
655 * then, in case of minimization, the corresponding row in the matrix
656 * will be
657 * a c + a e(y)
658 * with a d = m, the (updated) common denominator of the matrix.
659 * In case of maximization, the row will be
660 * -a c - a e(y)
662 static void sol_add(struct isl_sol *sol, struct isl_tab *tab)
664 struct isl_basic_set *bset = NULL;
665 struct isl_mat *mat = NULL;
666 unsigned off;
667 int row;
668 isl_int m;
670 if (sol->error || !tab)
671 goto error;
673 if (tab->empty && !sol->add_empty)
674 return;
675 if (sol->context->op->is_empty(sol->context))
676 return;
678 bset = sol_domain(sol);
680 if (tab->empty) {
681 sol_push_sol(sol, bset, NULL);
682 return;
685 off = 2 + tab->M;
687 mat = isl_mat_alloc(tab->mat->ctx, 1 + sol->n_out,
688 1 + tab->n_param + tab->n_div);
689 if (!mat)
690 goto error;
692 isl_int_init(m);
694 isl_seq_clr(mat->row[0] + 1, mat->n_col - 1);
695 isl_int_set_si(mat->row[0][0], 1);
696 for (row = 0; row < sol->n_out; ++row) {
697 int i = tab->n_param + row;
698 int r, j;
700 isl_seq_clr(mat->row[1 + row], mat->n_col);
701 if (!tab->var[i].is_row) {
702 if (tab->M)
703 isl_die(mat->ctx, isl_error_invalid,
704 "unbounded optimum", goto error2);
705 continue;
708 r = tab->var[i].index;
709 if (tab->M &&
710 isl_int_ne(tab->mat->row[r][2], tab->mat->row[r][0]))
711 isl_die(mat->ctx, isl_error_invalid,
712 "unbounded optimum", goto error2);
713 isl_int_gcd(m, mat->row[0][0], tab->mat->row[r][0]);
714 isl_int_divexact(m, tab->mat->row[r][0], m);
715 scale_rows(mat, m, 1 + row);
716 isl_int_divexact(m, mat->row[0][0], tab->mat->row[r][0]);
717 isl_int_mul(mat->row[1 + row][0], m, tab->mat->row[r][1]);
718 for (j = 0; j < tab->n_param; ++j) {
719 int col;
720 if (tab->var[j].is_row)
721 continue;
722 col = tab->var[j].index;
723 isl_int_mul(mat->row[1 + row][1 + j], m,
724 tab->mat->row[r][off + col]);
726 for (j = 0; j < tab->n_div; ++j) {
727 int col;
728 if (tab->var[tab->n_var - tab->n_div+j].is_row)
729 continue;
730 col = tab->var[tab->n_var - tab->n_div+j].index;
731 isl_int_mul(mat->row[1 + row][1 + tab->n_param + j], m,
732 tab->mat->row[r][off + col]);
734 if (sol->max)
735 isl_seq_neg(mat->row[1 + row], mat->row[1 + row],
736 mat->n_col);
739 isl_int_clear(m);
741 sol_push_sol_mat(sol, bset, mat);
742 return;
743 error2:
744 isl_int_clear(m);
745 error:
746 isl_basic_set_free(bset);
747 isl_mat_free(mat);
748 sol->error = 1;
751 struct isl_sol_map {
752 struct isl_sol sol;
753 struct isl_map *map;
754 struct isl_set *empty;
757 static void sol_map_free(struct isl_sol *sol)
759 struct isl_sol_map *sol_map = (struct isl_sol_map *) sol;
760 isl_map_free(sol_map->map);
761 isl_set_free(sol_map->empty);
764 /* This function is called for parts of the context where there is
765 * no solution, with "bset" corresponding to the context tableau.
766 * Simply add the basic set to the set "empty".
768 static void sol_map_add_empty(struct isl_sol_map *sol,
769 struct isl_basic_set *bset)
771 if (!bset || !sol->empty)
772 goto error;
774 sol->empty = isl_set_grow(sol->empty, 1);
775 bset = isl_basic_set_simplify(bset);
776 bset = isl_basic_set_finalize(bset);
777 sol->empty = isl_set_add_basic_set(sol->empty, isl_basic_set_copy(bset));
778 if (!sol->empty)
779 goto error;
780 isl_basic_set_free(bset);
781 return;
782 error:
783 isl_basic_set_free(bset);
784 sol->sol.error = 1;
787 static void sol_map_add_empty_wrap(struct isl_sol *sol,
788 struct isl_basic_set *bset)
790 sol_map_add_empty((struct isl_sol_map *)sol, bset);
793 /* Given a basic set "dom" that represents the context and a tuple of
794 * affine expressions "ma" defined over this domain, construct a basic map
795 * that expresses this function on the domain.
797 static void sol_map_add(struct isl_sol_map *sol,
798 __isl_take isl_basic_set *dom, __isl_take isl_multi_aff *ma)
800 isl_basic_map *bmap;
802 if (sol->sol.error || !dom || !ma)
803 goto error;
805 bmap = isl_basic_map_from_multi_aff2(ma, sol->sol.rational);
806 bmap = isl_basic_map_intersect_domain(bmap, dom);
807 sol->map = isl_map_grow(sol->map, 1);
808 sol->map = isl_map_add_basic_map(sol->map, bmap);
809 if (!sol->map)
810 sol->sol.error = 1;
811 return;
812 error:
813 isl_basic_set_free(dom);
814 isl_multi_aff_free(ma);
815 sol->sol.error = 1;
818 static void sol_map_add_wrap(struct isl_sol *sol,
819 __isl_take isl_basic_set *dom, __isl_take isl_multi_aff *ma)
821 sol_map_add((struct isl_sol_map *)sol, dom, ma);
825 /* Store the "parametric constant" of row "row" of tableau "tab" in "line",
826 * i.e., the constant term and the coefficients of all variables that
827 * appear in the context tableau.
828 * Note that the coefficient of the big parameter M is NOT copied.
829 * The context tableau may not have a big parameter and even when it
830 * does, it is a different big parameter.
832 static void get_row_parameter_line(struct isl_tab *tab, int row, isl_int *line)
834 int i;
835 unsigned off = 2 + tab->M;
837 isl_int_set(line[0], tab->mat->row[row][1]);
838 for (i = 0; i < tab->n_param; ++i) {
839 if (tab->var[i].is_row)
840 isl_int_set_si(line[1 + i], 0);
841 else {
842 int col = tab->var[i].index;
843 isl_int_set(line[1 + i], tab->mat->row[row][off + col]);
846 for (i = 0; i < tab->n_div; ++i) {
847 if (tab->var[tab->n_var - tab->n_div + i].is_row)
848 isl_int_set_si(line[1 + tab->n_param + i], 0);
849 else {
850 int col = tab->var[tab->n_var - tab->n_div + i].index;
851 isl_int_set(line[1 + tab->n_param + i],
852 tab->mat->row[row][off + col]);
857 /* Check if rows "row1" and "row2" have identical "parametric constants",
858 * as explained above.
859 * In this case, we also insist that the coefficients of the big parameter
860 * be the same as the values of the constants will only be the same
861 * if these coefficients are also the same.
863 static int identical_parameter_line(struct isl_tab *tab, int row1, int row2)
865 int i;
866 unsigned off = 2 + tab->M;
868 if (isl_int_ne(tab->mat->row[row1][1], tab->mat->row[row2][1]))
869 return 0;
871 if (tab->M && isl_int_ne(tab->mat->row[row1][2],
872 tab->mat->row[row2][2]))
873 return 0;
875 for (i = 0; i < tab->n_param + tab->n_div; ++i) {
876 int pos = i < tab->n_param ? i :
877 tab->n_var - tab->n_div + i - tab->n_param;
878 int col;
880 if (tab->var[pos].is_row)
881 continue;
882 col = tab->var[pos].index;
883 if (isl_int_ne(tab->mat->row[row1][off + col],
884 tab->mat->row[row2][off + col]))
885 return 0;
887 return 1;
890 /* Return an inequality that expresses that the "parametric constant"
891 * should be non-negative.
892 * This function is only called when the coefficient of the big parameter
893 * is equal to zero.
895 static struct isl_vec *get_row_parameter_ineq(struct isl_tab *tab, int row)
897 struct isl_vec *ineq;
899 ineq = isl_vec_alloc(tab->mat->ctx, 1 + tab->n_param + tab->n_div);
900 if (!ineq)
901 return NULL;
903 get_row_parameter_line(tab, row, ineq->el);
904 if (ineq)
905 ineq = isl_vec_normalize(ineq);
907 return ineq;
910 /* Normalize a div expression of the form
912 * [(g*f(x) + c)/(g * m)]
914 * with c the constant term and f(x) the remaining coefficients, to
916 * [(f(x) + [c/g])/m]
918 static void normalize_div(__isl_keep isl_vec *div)
920 isl_ctx *ctx = isl_vec_get_ctx(div);
921 int len = div->size - 2;
923 isl_seq_gcd(div->el + 2, len, &ctx->normalize_gcd);
924 isl_int_gcd(ctx->normalize_gcd, ctx->normalize_gcd, div->el[0]);
926 if (isl_int_is_one(ctx->normalize_gcd))
927 return;
929 isl_int_divexact(div->el[0], div->el[0], ctx->normalize_gcd);
930 isl_int_fdiv_q(div->el[1], div->el[1], ctx->normalize_gcd);
931 isl_seq_scale_down(div->el + 2, div->el + 2, ctx->normalize_gcd, len);
934 /* Return an integer division for use in a parametric cut based
935 * on the given row.
936 * In particular, let the parametric constant of the row be
938 * \sum_i a_i y_i
940 * where y_0 = 1, but none of the y_i corresponds to the big parameter M.
941 * The div returned is equal to
943 * floor(\sum_i {-a_i} y_i) = floor((\sum_i (-a_i mod d) y_i)/d)
945 static struct isl_vec *get_row_parameter_div(struct isl_tab *tab, int row)
947 struct isl_vec *div;
949 div = isl_vec_alloc(tab->mat->ctx, 1 + 1 + tab->n_param + tab->n_div);
950 if (!div)
951 return NULL;
953 isl_int_set(div->el[0], tab->mat->row[row][0]);
954 get_row_parameter_line(tab, row, div->el + 1);
955 isl_seq_neg(div->el + 1, div->el + 1, div->size - 1);
956 normalize_div(div);
957 isl_seq_fdiv_r(div->el + 1, div->el + 1, div->el[0], div->size - 1);
959 return div;
962 /* Return an integer division for use in transferring an integrality constraint
963 * to the context.
964 * In particular, let the parametric constant of the row be
966 * \sum_i a_i y_i
968 * where y_0 = 1, but none of the y_i corresponds to the big parameter M.
969 * The the returned div is equal to
971 * floor(\sum_i {a_i} y_i) = floor((\sum_i (a_i mod d) y_i)/d)
973 static struct isl_vec *get_row_split_div(struct isl_tab *tab, int row)
975 struct isl_vec *div;
977 div = isl_vec_alloc(tab->mat->ctx, 1 + 1 + tab->n_param + tab->n_div);
978 if (!div)
979 return NULL;
981 isl_int_set(div->el[0], tab->mat->row[row][0]);
982 get_row_parameter_line(tab, row, div->el + 1);
983 normalize_div(div);
984 isl_seq_fdiv_r(div->el + 1, div->el + 1, div->el[0], div->size - 1);
986 return div;
989 /* Construct and return an inequality that expresses an upper bound
990 * on the given div.
991 * In particular, if the div is given by
993 * d = floor(e/m)
995 * then the inequality expresses
997 * m d <= e
999 static __isl_give isl_vec *ineq_for_div(__isl_keep isl_basic_set *bset,
1000 unsigned div)
1002 unsigned total;
1003 unsigned div_pos;
1004 struct isl_vec *ineq;
1006 if (!bset)
1007 return NULL;
1009 total = isl_basic_set_total_dim(bset);
1010 div_pos = 1 + total - bset->n_div + div;
1012 ineq = isl_vec_alloc(bset->ctx, 1 + total);
1013 if (!ineq)
1014 return NULL;
1016 isl_seq_cpy(ineq->el, bset->div[div] + 1, 1 + total);
1017 isl_int_neg(ineq->el[div_pos], bset->div[div][0]);
1018 return ineq;
1021 /* Given a row in the tableau and a div that was created
1022 * using get_row_split_div and that has been constrained to equality, i.e.,
1024 * d = floor(\sum_i {a_i} y_i) = \sum_i {a_i} y_i
1026 * replace the expression "\sum_i {a_i} y_i" in the row by d,
1027 * i.e., we subtract "\sum_i {a_i} y_i" and add 1 d.
1028 * The coefficients of the non-parameters in the tableau have been
1029 * verified to be integral. We can therefore simply replace coefficient b
1030 * by floor(b). For the coefficients of the parameters we have
1031 * floor(a_i) = a_i - {a_i}, while for the other coefficients, we have
1032 * floor(b) = b.
1034 static struct isl_tab *set_row_cst_to_div(struct isl_tab *tab, int row, int div)
1036 isl_seq_fdiv_q(tab->mat->row[row] + 1, tab->mat->row[row] + 1,
1037 tab->mat->row[row][0], 1 + tab->M + tab->n_col);
1039 isl_int_set_si(tab->mat->row[row][0], 1);
1041 if (tab->var[tab->n_var - tab->n_div + div].is_row) {
1042 int drow = tab->var[tab->n_var - tab->n_div + div].index;
1044 isl_assert(tab->mat->ctx,
1045 isl_int_is_one(tab->mat->row[drow][0]), goto error);
1046 isl_seq_combine(tab->mat->row[row] + 1,
1047 tab->mat->ctx->one, tab->mat->row[row] + 1,
1048 tab->mat->ctx->one, tab->mat->row[drow] + 1,
1049 1 + tab->M + tab->n_col);
1050 } else {
1051 int dcol = tab->var[tab->n_var - tab->n_div + div].index;
1053 isl_int_add_ui(tab->mat->row[row][2 + tab->M + dcol],
1054 tab->mat->row[row][2 + tab->M + dcol], 1);
1057 return tab;
1058 error:
1059 isl_tab_free(tab);
1060 return NULL;
1063 /* Check if the (parametric) constant of the given row is obviously
1064 * negative, meaning that we don't need to consult the context tableau.
1065 * If there is a big parameter and its coefficient is non-zero,
1066 * then this coefficient determines the outcome.
1067 * Otherwise, we check whether the constant is negative and
1068 * all non-zero coefficients of parameters are negative and
1069 * belong to non-negative parameters.
1071 static int is_obviously_neg(struct isl_tab *tab, int row)
1073 int i;
1074 int col;
1075 unsigned off = 2 + tab->M;
1077 if (tab->M) {
1078 if (isl_int_is_pos(tab->mat->row[row][2]))
1079 return 0;
1080 if (isl_int_is_neg(tab->mat->row[row][2]))
1081 return 1;
1084 if (isl_int_is_nonneg(tab->mat->row[row][1]))
1085 return 0;
1086 for (i = 0; i < tab->n_param; ++i) {
1087 /* Eliminated parameter */
1088 if (tab->var[i].is_row)
1089 continue;
1090 col = tab->var[i].index;
1091 if (isl_int_is_zero(tab->mat->row[row][off + col]))
1092 continue;
1093 if (!tab->var[i].is_nonneg)
1094 return 0;
1095 if (isl_int_is_pos(tab->mat->row[row][off + col]))
1096 return 0;
1098 for (i = 0; i < tab->n_div; ++i) {
1099 if (tab->var[tab->n_var - tab->n_div + i].is_row)
1100 continue;
1101 col = tab->var[tab->n_var - tab->n_div + i].index;
1102 if (isl_int_is_zero(tab->mat->row[row][off + col]))
1103 continue;
1104 if (!tab->var[tab->n_var - tab->n_div + i].is_nonneg)
1105 return 0;
1106 if (isl_int_is_pos(tab->mat->row[row][off + col]))
1107 return 0;
1109 return 1;
1112 /* Check if the (parametric) constant of the given row is obviously
1113 * non-negative, meaning that we don't need to consult the context tableau.
1114 * If there is a big parameter and its coefficient is non-zero,
1115 * then this coefficient determines the outcome.
1116 * Otherwise, we check whether the constant is non-negative and
1117 * all non-zero coefficients of parameters are positive and
1118 * belong to non-negative parameters.
1120 static int is_obviously_nonneg(struct isl_tab *tab, int row)
1122 int i;
1123 int col;
1124 unsigned off = 2 + tab->M;
1126 if (tab->M) {
1127 if (isl_int_is_pos(tab->mat->row[row][2]))
1128 return 1;
1129 if (isl_int_is_neg(tab->mat->row[row][2]))
1130 return 0;
1133 if (isl_int_is_neg(tab->mat->row[row][1]))
1134 return 0;
1135 for (i = 0; i < tab->n_param; ++i) {
1136 /* Eliminated parameter */
1137 if (tab->var[i].is_row)
1138 continue;
1139 col = tab->var[i].index;
1140 if (isl_int_is_zero(tab->mat->row[row][off + col]))
1141 continue;
1142 if (!tab->var[i].is_nonneg)
1143 return 0;
1144 if (isl_int_is_neg(tab->mat->row[row][off + col]))
1145 return 0;
1147 for (i = 0; i < tab->n_div; ++i) {
1148 if (tab->var[tab->n_var - tab->n_div + i].is_row)
1149 continue;
1150 col = tab->var[tab->n_var - tab->n_div + i].index;
1151 if (isl_int_is_zero(tab->mat->row[row][off + col]))
1152 continue;
1153 if (!tab->var[tab->n_var - tab->n_div + i].is_nonneg)
1154 return 0;
1155 if (isl_int_is_neg(tab->mat->row[row][off + col]))
1156 return 0;
1158 return 1;
1161 /* Given a row r and two columns, return the column that would
1162 * lead to the lexicographically smallest increment in the sample
1163 * solution when leaving the basis in favor of the row.
1164 * Pivoting with column c will increment the sample value by a non-negative
1165 * constant times a_{V,c}/a_{r,c}, with a_{V,c} the elements of column c
1166 * corresponding to the non-parametric variables.
1167 * If variable v appears in a column c_v, then a_{v,c} = 1 iff c = c_v,
1168 * with all other entries in this virtual row equal to zero.
1169 * If variable v appears in a row, then a_{v,c} is the element in column c
1170 * of that row.
1172 * Let v be the first variable with a_{v,c1}/a_{r,c1} != a_{v,c2}/a_{r,c2}.
1173 * Then if a_{v,c1}/a_{r,c1} < a_{v,c2}/a_{r,c2}, i.e.,
1174 * a_{v,c2} a_{r,c1} - a_{v,c1} a_{r,c2} > 0, c1 results in the minimal
1175 * increment. Otherwise, it's c2.
1177 static int lexmin_col_pair(struct isl_tab *tab,
1178 int row, int col1, int col2, isl_int tmp)
1180 int i;
1181 isl_int *tr;
1183 tr = tab->mat->row[row] + 2 + tab->M;
1185 for (i = tab->n_param; i < tab->n_var - tab->n_div; ++i) {
1186 int s1, s2;
1187 isl_int *r;
1189 if (!tab->var[i].is_row) {
1190 if (tab->var[i].index == col1)
1191 return col2;
1192 if (tab->var[i].index == col2)
1193 return col1;
1194 continue;
1197 if (tab->var[i].index == row)
1198 continue;
1200 r = tab->mat->row[tab->var[i].index] + 2 + tab->M;
1201 s1 = isl_int_sgn(r[col1]);
1202 s2 = isl_int_sgn(r[col2]);
1203 if (s1 == 0 && s2 == 0)
1204 continue;
1205 if (s1 < s2)
1206 return col1;
1207 if (s2 < s1)
1208 return col2;
1210 isl_int_mul(tmp, r[col2], tr[col1]);
1211 isl_int_submul(tmp, r[col1], tr[col2]);
1212 if (isl_int_is_pos(tmp))
1213 return col1;
1214 if (isl_int_is_neg(tmp))
1215 return col2;
1217 return -1;
1220 /* Does the index into the tab->var or tab->con array "index"
1221 * correspond to a variable in the context tableau?
1222 * In particular, it needs to be an index into the tab->var array and
1223 * it needs to refer to either one of the first tab->n_param variables or
1224 * one of the last tab->n_div variables.
1226 static int is_parameter_var(struct isl_tab *tab, int index)
1228 if (index < 0)
1229 return 0;
1230 if (index < tab->n_param)
1231 return 1;
1232 if (index >= tab->n_var - tab->n_div)
1233 return 1;
1234 return 0;
1237 /* Does column "col" of "tab" refer to a variable in the context tableau?
1239 static int col_is_parameter_var(struct isl_tab *tab, int col)
1241 return is_parameter_var(tab, tab->col_var[col]);
1244 /* Does row "row" of "tab" refer to a variable in the context tableau?
1246 static int row_is_parameter_var(struct isl_tab *tab, int row)
1248 return is_parameter_var(tab, tab->row_var[row]);
1251 /* Given a row in the tableau, find and return the column that would
1252 * result in the lexicographically smallest, but positive, increment
1253 * in the sample point.
1254 * If there is no such column, then return tab->n_col.
1255 * If anything goes wrong, return -1.
1257 static int lexmin_pivot_col(struct isl_tab *tab, int row)
1259 int j;
1260 int col = tab->n_col;
1261 isl_int *tr;
1262 isl_int tmp;
1264 tr = tab->mat->row[row] + 2 + tab->M;
1266 isl_int_init(tmp);
1268 for (j = tab->n_dead; j < tab->n_col; ++j) {
1269 if (col_is_parameter_var(tab, j))
1270 continue;
1272 if (!isl_int_is_pos(tr[j]))
1273 continue;
1275 if (col == tab->n_col)
1276 col = j;
1277 else
1278 col = lexmin_col_pair(tab, row, col, j, tmp);
1279 isl_assert(tab->mat->ctx, col >= 0, goto error);
1282 isl_int_clear(tmp);
1283 return col;
1284 error:
1285 isl_int_clear(tmp);
1286 return -1;
1289 /* Return the first known violated constraint, i.e., a non-negative
1290 * constraint that currently has an either obviously negative value
1291 * or a previously determined to be negative value.
1293 * If any constraint has a negative coefficient for the big parameter,
1294 * if any, then we return one of these first.
1296 static int first_neg(struct isl_tab *tab)
1298 int row;
1300 if (tab->M)
1301 for (row = tab->n_redundant; row < tab->n_row; ++row) {
1302 if (!isl_tab_var_from_row(tab, row)->is_nonneg)
1303 continue;
1304 if (!isl_int_is_neg(tab->mat->row[row][2]))
1305 continue;
1306 if (tab->row_sign)
1307 tab->row_sign[row] = isl_tab_row_neg;
1308 return row;
1310 for (row = tab->n_redundant; row < tab->n_row; ++row) {
1311 if (!isl_tab_var_from_row(tab, row)->is_nonneg)
1312 continue;
1313 if (tab->row_sign) {
1314 if (tab->row_sign[row] == 0 &&
1315 is_obviously_neg(tab, row))
1316 tab->row_sign[row] = isl_tab_row_neg;
1317 if (tab->row_sign[row] != isl_tab_row_neg)
1318 continue;
1319 } else if (!is_obviously_neg(tab, row))
1320 continue;
1321 return row;
1323 return -1;
1326 /* Check whether the invariant that all columns are lexico-positive
1327 * is satisfied. This function is not called from the current code
1328 * but is useful during debugging.
1330 static void check_lexpos(struct isl_tab *tab) __attribute__ ((unused));
1331 static void check_lexpos(struct isl_tab *tab)
1333 unsigned off = 2 + tab->M;
1334 int col;
1335 int var;
1336 int row;
1338 for (col = tab->n_dead; col < tab->n_col; ++col) {
1339 if (col_is_parameter_var(tab, col))
1340 continue;
1341 for (var = tab->n_param; var < tab->n_var - tab->n_div; ++var) {
1342 if (!tab->var[var].is_row) {
1343 if (tab->var[var].index == col)
1344 break;
1345 else
1346 continue;
1348 row = tab->var[var].index;
1349 if (isl_int_is_zero(tab->mat->row[row][off + col]))
1350 continue;
1351 if (isl_int_is_pos(tab->mat->row[row][off + col]))
1352 break;
1353 fprintf(stderr, "lexneg column %d (row %d)\n",
1354 col, row);
1356 if (var >= tab->n_var - tab->n_div)
1357 fprintf(stderr, "zero column %d\n", col);
1361 /* Report to the caller that the given constraint is part of an encountered
1362 * conflict.
1364 static int report_conflicting_constraint(struct isl_tab *tab, int con)
1366 return tab->conflict(con, tab->conflict_user);
1369 /* Given a conflicting row in the tableau, report all constraints
1370 * involved in the row to the caller. That is, the row itself
1371 * (if it represents a constraint) and all constraint columns with
1372 * non-zero (and therefore negative) coefficients.
1374 static int report_conflict(struct isl_tab *tab, int row)
1376 int j;
1377 isl_int *tr;
1379 if (!tab->conflict)
1380 return 0;
1382 if (tab->row_var[row] < 0 &&
1383 report_conflicting_constraint(tab, ~tab->row_var[row]) < 0)
1384 return -1;
1386 tr = tab->mat->row[row] + 2 + tab->M;
1388 for (j = tab->n_dead; j < tab->n_col; ++j) {
1389 if (col_is_parameter_var(tab, j))
1390 continue;
1392 if (!isl_int_is_neg(tr[j]))
1393 continue;
1395 if (tab->col_var[j] < 0 &&
1396 report_conflicting_constraint(tab, ~tab->col_var[j]) < 0)
1397 return -1;
1400 return 0;
1403 /* Resolve all known or obviously violated constraints through pivoting.
1404 * In particular, as long as we can find any violated constraint, we
1405 * look for a pivoting column that would result in the lexicographically
1406 * smallest increment in the sample point. If there is no such column
1407 * then the tableau is infeasible.
1409 static int restore_lexmin(struct isl_tab *tab) WARN_UNUSED;
1410 static int restore_lexmin(struct isl_tab *tab)
1412 int row, col;
1414 if (!tab)
1415 return -1;
1416 if (tab->empty)
1417 return 0;
1418 while ((row = first_neg(tab)) != -1) {
1419 col = lexmin_pivot_col(tab, row);
1420 if (col >= tab->n_col) {
1421 if (report_conflict(tab, row) < 0)
1422 return -1;
1423 if (isl_tab_mark_empty(tab) < 0)
1424 return -1;
1425 return 0;
1427 if (col < 0)
1428 return -1;
1429 if (isl_tab_pivot(tab, row, col) < 0)
1430 return -1;
1432 return 0;
1435 /* Given a row that represents an equality, look for an appropriate
1436 * pivoting column.
1437 * In particular, if there are any non-zero coefficients among
1438 * the non-parameter variables, then we take the last of these
1439 * variables. Eliminating this variable in terms of the other
1440 * variables and/or parameters does not influence the property
1441 * that all column in the initial tableau are lexicographically
1442 * positive. The row corresponding to the eliminated variable
1443 * will only have non-zero entries below the diagonal of the
1444 * initial tableau. That is, we transform
1446 * I I
1447 * 1 into a
1448 * I I
1450 * If there is no such non-parameter variable, then we are dealing with
1451 * pure parameter equality and we pick any parameter with coefficient 1 or -1
1452 * for elimination. This will ensure that the eliminated parameter
1453 * always has an integer value whenever all the other parameters are integral.
1454 * If there is no such parameter then we return -1.
1456 static int last_var_col_or_int_par_col(struct isl_tab *tab, int row)
1458 unsigned off = 2 + tab->M;
1459 int i;
1461 for (i = tab->n_var - tab->n_div - 1; i >= 0 && i >= tab->n_param; --i) {
1462 int col;
1463 if (tab->var[i].is_row)
1464 continue;
1465 col = tab->var[i].index;
1466 if (col <= tab->n_dead)
1467 continue;
1468 if (!isl_int_is_zero(tab->mat->row[row][off + col]))
1469 return col;
1471 for (i = tab->n_dead; i < tab->n_col; ++i) {
1472 if (isl_int_is_one(tab->mat->row[row][off + i]))
1473 return i;
1474 if (isl_int_is_negone(tab->mat->row[row][off + i]))
1475 return i;
1477 return -1;
1480 /* Add an equality that is known to be valid to the tableau.
1481 * We first check if we can eliminate a variable or a parameter.
1482 * If not, we add the equality as two inequalities.
1483 * In this case, the equality was a pure parameter equality and there
1484 * is no need to resolve any constraint violations.
1486 * This function assumes that at least two more rows and at least
1487 * two more elements in the constraint array are available in the tableau.
1489 static struct isl_tab *add_lexmin_valid_eq(struct isl_tab *tab, isl_int *eq)
1491 int i;
1492 int r;
1494 if (!tab)
1495 return NULL;
1496 r = isl_tab_add_row(tab, eq);
1497 if (r < 0)
1498 goto error;
1500 r = tab->con[r].index;
1501 i = last_var_col_or_int_par_col(tab, r);
1502 if (i < 0) {
1503 tab->con[r].is_nonneg = 1;
1504 if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r]) < 0)
1505 goto error;
1506 isl_seq_neg(eq, eq, 1 + tab->n_var);
1507 r = isl_tab_add_row(tab, eq);
1508 if (r < 0)
1509 goto error;
1510 tab->con[r].is_nonneg = 1;
1511 if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r]) < 0)
1512 goto error;
1513 } else {
1514 if (isl_tab_pivot(tab, r, i) < 0)
1515 goto error;
1516 if (isl_tab_kill_col(tab, i) < 0)
1517 goto error;
1518 tab->n_eq++;
1521 return tab;
1522 error:
1523 isl_tab_free(tab);
1524 return NULL;
1527 /* Check if the given row is a pure constant.
1529 static int is_constant(struct isl_tab *tab, int row)
1531 unsigned off = 2 + tab->M;
1533 return isl_seq_first_non_zero(tab->mat->row[row] + off + tab->n_dead,
1534 tab->n_col - tab->n_dead) == -1;
1537 /* Is the given row a parametric constant?
1538 * That is, does it only involve variables that also appear in the context?
1540 static int is_parametric_constant(struct isl_tab *tab, int row)
1542 unsigned off = 2 + tab->M;
1543 int col;
1545 for (col = tab->n_dead; col < tab->n_col; ++col) {
1546 if (col_is_parameter_var(tab, col))
1547 continue;
1548 if (isl_int_is_zero(tab->mat->row[row][off + col]))
1549 continue;
1550 return 0;
1553 return 1;
1556 /* Add an equality that may or may not be valid to the tableau.
1557 * If the resulting row is a pure constant, then it must be zero.
1558 * Otherwise, the resulting tableau is empty.
1560 * If the row is not a pure constant, then we add two inequalities,
1561 * each time checking that they can be satisfied.
1562 * In the end we try to use one of the two constraints to eliminate
1563 * a column.
1565 * This function assumes that at least two more rows and at least
1566 * two more elements in the constraint array are available in the tableau.
1568 static int add_lexmin_eq(struct isl_tab *tab, isl_int *eq) WARN_UNUSED;
1569 static int add_lexmin_eq(struct isl_tab *tab, isl_int *eq)
1571 int r1, r2;
1572 int row;
1573 struct isl_tab_undo *snap;
1575 if (!tab)
1576 return -1;
1577 snap = isl_tab_snap(tab);
1578 r1 = isl_tab_add_row(tab, eq);
1579 if (r1 < 0)
1580 return -1;
1581 tab->con[r1].is_nonneg = 1;
1582 if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r1]) < 0)
1583 return -1;
1585 row = tab->con[r1].index;
1586 if (is_constant(tab, row)) {
1587 if (!isl_int_is_zero(tab->mat->row[row][1]) ||
1588 (tab->M && !isl_int_is_zero(tab->mat->row[row][2]))) {
1589 if (isl_tab_mark_empty(tab) < 0)
1590 return -1;
1591 return 0;
1593 if (isl_tab_rollback(tab, snap) < 0)
1594 return -1;
1595 return 0;
1598 if (restore_lexmin(tab) < 0)
1599 return -1;
1600 if (tab->empty)
1601 return 0;
1603 isl_seq_neg(eq, eq, 1 + tab->n_var);
1605 r2 = isl_tab_add_row(tab, eq);
1606 if (r2 < 0)
1607 return -1;
1608 tab->con[r2].is_nonneg = 1;
1609 if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r2]) < 0)
1610 return -1;
1612 if (restore_lexmin(tab) < 0)
1613 return -1;
1614 if (tab->empty)
1615 return 0;
1617 if (!tab->con[r1].is_row) {
1618 if (isl_tab_kill_col(tab, tab->con[r1].index) < 0)
1619 return -1;
1620 } else if (!tab->con[r2].is_row) {
1621 if (isl_tab_kill_col(tab, tab->con[r2].index) < 0)
1622 return -1;
1625 if (tab->bmap) {
1626 tab->bmap = isl_basic_map_add_ineq(tab->bmap, eq);
1627 if (isl_tab_push(tab, isl_tab_undo_bmap_ineq) < 0)
1628 return -1;
1629 isl_seq_neg(eq, eq, 1 + tab->n_var);
1630 tab->bmap = isl_basic_map_add_ineq(tab->bmap, eq);
1631 isl_seq_neg(eq, eq, 1 + tab->n_var);
1632 if (isl_tab_push(tab, isl_tab_undo_bmap_ineq) < 0)
1633 return -1;
1634 if (!tab->bmap)
1635 return -1;
1638 return 0;
1641 /* Add an inequality to the tableau, resolving violations using
1642 * restore_lexmin.
1644 * This function assumes that at least one more row and at least
1645 * one more element in the constraint array are available in the tableau.
1647 static struct isl_tab *add_lexmin_ineq(struct isl_tab *tab, isl_int *ineq)
1649 int r;
1651 if (!tab)
1652 return NULL;
1653 if (tab->bmap) {
1654 tab->bmap = isl_basic_map_add_ineq(tab->bmap, ineq);
1655 if (isl_tab_push(tab, isl_tab_undo_bmap_ineq) < 0)
1656 goto error;
1657 if (!tab->bmap)
1658 goto error;
1660 r = isl_tab_add_row(tab, ineq);
1661 if (r < 0)
1662 goto error;
1663 tab->con[r].is_nonneg = 1;
1664 if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r]) < 0)
1665 goto error;
1666 if (isl_tab_row_is_redundant(tab, tab->con[r].index)) {
1667 if (isl_tab_mark_redundant(tab, tab->con[r].index) < 0)
1668 goto error;
1669 return tab;
1672 if (restore_lexmin(tab) < 0)
1673 goto error;
1674 if (!tab->empty && tab->con[r].is_row &&
1675 isl_tab_row_is_redundant(tab, tab->con[r].index))
1676 if (isl_tab_mark_redundant(tab, tab->con[r].index) < 0)
1677 goto error;
1678 return tab;
1679 error:
1680 isl_tab_free(tab);
1681 return NULL;
1684 /* Check if the coefficients of the parameters are all integral.
1686 static int integer_parameter(struct isl_tab *tab, int row)
1688 int i;
1689 int col;
1690 unsigned off = 2 + tab->M;
1692 for (i = 0; i < tab->n_param; ++i) {
1693 /* Eliminated parameter */
1694 if (tab->var[i].is_row)
1695 continue;
1696 col = tab->var[i].index;
1697 if (!isl_int_is_divisible_by(tab->mat->row[row][off + col],
1698 tab->mat->row[row][0]))
1699 return 0;
1701 for (i = 0; i < tab->n_div; ++i) {
1702 if (tab->var[tab->n_var - tab->n_div + i].is_row)
1703 continue;
1704 col = tab->var[tab->n_var - tab->n_div + i].index;
1705 if (!isl_int_is_divisible_by(tab->mat->row[row][off + col],
1706 tab->mat->row[row][0]))
1707 return 0;
1709 return 1;
1712 /* Check if the coefficients of the non-parameter variables are all integral.
1714 static int integer_variable(struct isl_tab *tab, int row)
1716 int i;
1717 unsigned off = 2 + tab->M;
1719 for (i = tab->n_dead; i < tab->n_col; ++i) {
1720 if (col_is_parameter_var(tab, i))
1721 continue;
1722 if (!isl_int_is_divisible_by(tab->mat->row[row][off + i],
1723 tab->mat->row[row][0]))
1724 return 0;
1726 return 1;
1729 /* Check if the constant term is integral.
1731 static int integer_constant(struct isl_tab *tab, int row)
1733 return isl_int_is_divisible_by(tab->mat->row[row][1],
1734 tab->mat->row[row][0]);
1737 #define I_CST 1 << 0
1738 #define I_PAR 1 << 1
1739 #define I_VAR 1 << 2
1741 /* Check for next (non-parameter) variable after "var" (first if var == -1)
1742 * that is non-integer and therefore requires a cut and return
1743 * the index of the variable.
1744 * For parametric tableaus, there are three parts in a row,
1745 * the constant, the coefficients of the parameters and the rest.
1746 * For each part, we check whether the coefficients in that part
1747 * are all integral and if so, set the corresponding flag in *f.
1748 * If the constant and the parameter part are integral, then the
1749 * current sample value is integral and no cut is required
1750 * (irrespective of whether the variable part is integral).
1752 static int next_non_integer_var(struct isl_tab *tab, int var, int *f)
1754 var = var < 0 ? tab->n_param : var + 1;
1756 for (; var < tab->n_var - tab->n_div; ++var) {
1757 int flags = 0;
1758 int row;
1759 if (!tab->var[var].is_row)
1760 continue;
1761 row = tab->var[var].index;
1762 if (integer_constant(tab, row))
1763 ISL_FL_SET(flags, I_CST);
1764 if (integer_parameter(tab, row))
1765 ISL_FL_SET(flags, I_PAR);
1766 if (ISL_FL_ISSET(flags, I_CST) && ISL_FL_ISSET(flags, I_PAR))
1767 continue;
1768 if (integer_variable(tab, row))
1769 ISL_FL_SET(flags, I_VAR);
1770 *f = flags;
1771 return var;
1773 return -1;
1776 /* Check for first (non-parameter) variable that is non-integer and
1777 * therefore requires a cut and return the corresponding row.
1778 * For parametric tableaus, there are three parts in a row,
1779 * the constant, the coefficients of the parameters and the rest.
1780 * For each part, we check whether the coefficients in that part
1781 * are all integral and if so, set the corresponding flag in *f.
1782 * If the constant and the parameter part are integral, then the
1783 * current sample value is integral and no cut is required
1784 * (irrespective of whether the variable part is integral).
1786 static int first_non_integer_row(struct isl_tab *tab, int *f)
1788 int var = next_non_integer_var(tab, -1, f);
1790 return var < 0 ? -1 : tab->var[var].index;
1793 /* Add a (non-parametric) cut to cut away the non-integral sample
1794 * value of the given row.
1796 * If the row is given by
1798 * m r = f + \sum_i a_i y_i
1800 * then the cut is
1802 * c = - {-f/m} + \sum_i {a_i/m} y_i >= 0
1804 * The big parameter, if any, is ignored, since it is assumed to be big
1805 * enough to be divisible by any integer.
1806 * If the tableau is actually a parametric tableau, then this function
1807 * is only called when all coefficients of the parameters are integral.
1808 * The cut therefore has zero coefficients for the parameters.
1810 * The current value is known to be negative, so row_sign, if it
1811 * exists, is set accordingly.
1813 * Return the row of the cut or -1.
1815 static int add_cut(struct isl_tab *tab, int row)
1817 int i;
1818 int r;
1819 isl_int *r_row;
1820 unsigned off = 2 + tab->M;
1822 if (isl_tab_extend_cons(tab, 1) < 0)
1823 return -1;
1824 r = isl_tab_allocate_con(tab);
1825 if (r < 0)
1826 return -1;
1828 r_row = tab->mat->row[tab->con[r].index];
1829 isl_int_set(r_row[0], tab->mat->row[row][0]);
1830 isl_int_neg(r_row[1], tab->mat->row[row][1]);
1831 isl_int_fdiv_r(r_row[1], r_row[1], tab->mat->row[row][0]);
1832 isl_int_neg(r_row[1], r_row[1]);
1833 if (tab->M)
1834 isl_int_set_si(r_row[2], 0);
1835 for (i = 0; i < tab->n_col; ++i)
1836 isl_int_fdiv_r(r_row[off + i],
1837 tab->mat->row[row][off + i], tab->mat->row[row][0]);
1839 tab->con[r].is_nonneg = 1;
1840 if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r]) < 0)
1841 return -1;
1842 if (tab->row_sign)
1843 tab->row_sign[tab->con[r].index] = isl_tab_row_neg;
1845 return tab->con[r].index;
1848 #define CUT_ALL 1
1849 #define CUT_ONE 0
1851 /* Given a non-parametric tableau, add cuts until an integer
1852 * sample point is obtained or until the tableau is determined
1853 * to be integer infeasible.
1854 * As long as there is any non-integer value in the sample point,
1855 * we add appropriate cuts, if possible, for each of these
1856 * non-integer values and then resolve the violated
1857 * cut constraints using restore_lexmin.
1858 * If one of the corresponding rows is equal to an integral
1859 * combination of variables/constraints plus a non-integral constant,
1860 * then there is no way to obtain an integer point and we return
1861 * a tableau that is marked empty.
1862 * The parameter cutting_strategy controls the strategy used when adding cuts
1863 * to remove non-integer points. CUT_ALL adds all possible cuts
1864 * before continuing the search. CUT_ONE adds only one cut at a time.
1866 static struct isl_tab *cut_to_integer_lexmin(struct isl_tab *tab,
1867 int cutting_strategy)
1869 int var;
1870 int row;
1871 int flags;
1873 if (!tab)
1874 return NULL;
1875 if (tab->empty)
1876 return tab;
1878 while ((var = next_non_integer_var(tab, -1, &flags)) != -1) {
1879 do {
1880 if (ISL_FL_ISSET(flags, I_VAR)) {
1881 if (isl_tab_mark_empty(tab) < 0)
1882 goto error;
1883 return tab;
1885 row = tab->var[var].index;
1886 row = add_cut(tab, row);
1887 if (row < 0)
1888 goto error;
1889 if (cutting_strategy == CUT_ONE)
1890 break;
1891 } while ((var = next_non_integer_var(tab, var, &flags)) != -1);
1892 if (restore_lexmin(tab) < 0)
1893 goto error;
1894 if (tab->empty)
1895 break;
1897 return tab;
1898 error:
1899 isl_tab_free(tab);
1900 return NULL;
1903 /* Check whether all the currently active samples also satisfy the inequality
1904 * "ineq" (treated as an equality if eq is set).
1905 * Remove those samples that do not.
1907 static struct isl_tab *check_samples(struct isl_tab *tab, isl_int *ineq, int eq)
1909 int i;
1910 isl_int v;
1912 if (!tab)
1913 return NULL;
1915 isl_assert(tab->mat->ctx, tab->bmap, goto error);
1916 isl_assert(tab->mat->ctx, tab->samples, goto error);
1917 isl_assert(tab->mat->ctx, tab->samples->n_col == 1 + tab->n_var, goto error);
1919 isl_int_init(v);
1920 for (i = tab->n_outside; i < tab->n_sample; ++i) {
1921 int sgn;
1922 isl_seq_inner_product(ineq, tab->samples->row[i],
1923 1 + tab->n_var, &v);
1924 sgn = isl_int_sgn(v);
1925 if (eq ? (sgn == 0) : (sgn >= 0))
1926 continue;
1927 tab = isl_tab_drop_sample(tab, i);
1928 if (!tab)
1929 break;
1931 isl_int_clear(v);
1933 return tab;
1934 error:
1935 isl_tab_free(tab);
1936 return NULL;
1939 /* Check whether the sample value of the tableau is finite,
1940 * i.e., either the tableau does not use a big parameter, or
1941 * all values of the variables are equal to the big parameter plus
1942 * some constant. This constant is the actual sample value.
1944 static int sample_is_finite(struct isl_tab *tab)
1946 int i;
1948 if (!tab->M)
1949 return 1;
1951 for (i = 0; i < tab->n_var; ++i) {
1952 int row;
1953 if (!tab->var[i].is_row)
1954 return 0;
1955 row = tab->var[i].index;
1956 if (isl_int_ne(tab->mat->row[row][0], tab->mat->row[row][2]))
1957 return 0;
1959 return 1;
1962 /* Check if the context tableau of sol has any integer points.
1963 * Leave tab in empty state if no integer point can be found.
1964 * If an integer point can be found and if moreover it is finite,
1965 * then it is added to the list of sample values.
1967 * This function is only called when none of the currently active sample
1968 * values satisfies the most recently added constraint.
1970 static struct isl_tab *check_integer_feasible(struct isl_tab *tab)
1972 struct isl_tab_undo *snap;
1974 if (!tab)
1975 return NULL;
1977 snap = isl_tab_snap(tab);
1978 if (isl_tab_push_basis(tab) < 0)
1979 goto error;
1981 tab = cut_to_integer_lexmin(tab, CUT_ALL);
1982 if (!tab)
1983 goto error;
1985 if (!tab->empty && sample_is_finite(tab)) {
1986 struct isl_vec *sample;
1988 sample = isl_tab_get_sample_value(tab);
1990 if (isl_tab_add_sample(tab, sample) < 0)
1991 goto error;
1994 if (!tab->empty && isl_tab_rollback(tab, snap) < 0)
1995 goto error;
1997 return tab;
1998 error:
1999 isl_tab_free(tab);
2000 return NULL;
2003 /* Check if any of the currently active sample values satisfies
2004 * the inequality "ineq" (an equality if eq is set).
2006 static int tab_has_valid_sample(struct isl_tab *tab, isl_int *ineq, int eq)
2008 int i;
2009 isl_int v;
2011 if (!tab)
2012 return -1;
2014 isl_assert(tab->mat->ctx, tab->bmap, return -1);
2015 isl_assert(tab->mat->ctx, tab->samples, return -1);
2016 isl_assert(tab->mat->ctx, tab->samples->n_col == 1 + tab->n_var, return -1);
2018 isl_int_init(v);
2019 for (i = tab->n_outside; i < tab->n_sample; ++i) {
2020 int sgn;
2021 isl_seq_inner_product(ineq, tab->samples->row[i],
2022 1 + tab->n_var, &v);
2023 sgn = isl_int_sgn(v);
2024 if (eq ? (sgn == 0) : (sgn >= 0))
2025 break;
2027 isl_int_clear(v);
2029 return i < tab->n_sample;
2032 /* Insert a div specified by "div" to the tableau "tab" at position "pos" and
2033 * return isl_bool_true if the div is obviously non-negative.
2035 static isl_bool context_tab_insert_div(struct isl_tab *tab, int pos,
2036 __isl_keep isl_vec *div,
2037 isl_stat (*add_ineq)(void *user, isl_int *), void *user)
2039 int i;
2040 int r;
2041 struct isl_mat *samples;
2042 int nonneg;
2044 r = isl_tab_insert_div(tab, pos, div, add_ineq, user);
2045 if (r < 0)
2046 return isl_bool_error;
2047 nonneg = tab->var[r].is_nonneg;
2048 tab->var[r].frozen = 1;
2050 samples = isl_mat_extend(tab->samples,
2051 tab->n_sample, 1 + tab->n_var);
2052 tab->samples = samples;
2053 if (!samples)
2054 return isl_bool_error;
2055 for (i = tab->n_outside; i < samples->n_row; ++i) {
2056 isl_seq_inner_product(div->el + 1, samples->row[i],
2057 div->size - 1, &samples->row[i][samples->n_col - 1]);
2058 isl_int_fdiv_q(samples->row[i][samples->n_col - 1],
2059 samples->row[i][samples->n_col - 1], div->el[0]);
2061 tab->samples = isl_mat_move_cols(tab->samples, 1 + pos,
2062 1 + tab->n_var - 1, 1);
2063 if (!tab->samples)
2064 return isl_bool_error;
2066 return nonneg;
2069 /* Add a div specified by "div" to both the main tableau and
2070 * the context tableau. In case of the main tableau, we only
2071 * need to add an extra div. In the context tableau, we also
2072 * need to express the meaning of the div.
2073 * Return the index of the div or -1 if anything went wrong.
2075 * The new integer division is added before any unknown integer
2076 * divisions in the context to ensure that it does not get
2077 * equated to some linear combination involving unknown integer
2078 * divisions.
2080 static int add_div(struct isl_tab *tab, struct isl_context *context,
2081 __isl_keep isl_vec *div)
2083 int r;
2084 int pos;
2085 isl_bool nonneg;
2086 struct isl_tab *context_tab = context->op->peek_tab(context);
2088 if (!tab || !context_tab)
2089 goto error;
2091 pos = context_tab->n_var - context->n_unknown;
2092 if ((nonneg = context->op->insert_div(context, pos, div)) < 0)
2093 goto error;
2095 if (!context->op->is_ok(context))
2096 goto error;
2098 pos = tab->n_var - context->n_unknown;
2099 if (isl_tab_extend_vars(tab, 1) < 0)
2100 goto error;
2101 r = isl_tab_insert_var(tab, pos);
2102 if (r < 0)
2103 goto error;
2104 if (nonneg)
2105 tab->var[r].is_nonneg = 1;
2106 tab->var[r].frozen = 1;
2107 tab->n_div++;
2109 return tab->n_div - 1 - context->n_unknown;
2110 error:
2111 context->op->invalidate(context);
2112 return -1;
2115 static int find_div(struct isl_tab *tab, isl_int *div, isl_int denom)
2117 int i;
2118 unsigned total = isl_basic_map_total_dim(tab->bmap);
2120 for (i = 0; i < tab->bmap->n_div; ++i) {
2121 if (isl_int_ne(tab->bmap->div[i][0], denom))
2122 continue;
2123 if (!isl_seq_eq(tab->bmap->div[i] + 1, div, 1 + total))
2124 continue;
2125 return i;
2127 return -1;
2130 /* Return the index of a div that corresponds to "div".
2131 * We first check if we already have such a div and if not, we create one.
2133 static int get_div(struct isl_tab *tab, struct isl_context *context,
2134 struct isl_vec *div)
2136 int d;
2137 struct isl_tab *context_tab = context->op->peek_tab(context);
2139 if (!context_tab)
2140 return -1;
2142 d = find_div(context_tab, div->el + 1, div->el[0]);
2143 if (d != -1)
2144 return d;
2146 return add_div(tab, context, div);
2149 /* Add a parametric cut to cut away the non-integral sample value
2150 * of the given row.
2151 * Let a_i be the coefficients of the constant term and the parameters
2152 * and let b_i be the coefficients of the variables or constraints
2153 * in basis of the tableau.
2154 * Let q be the div q = floor(\sum_i {-a_i} y_i).
2156 * The cut is expressed as
2158 * c = \sum_i -{-a_i} y_i + \sum_i {b_i} x_i + q >= 0
2160 * If q did not already exist in the context tableau, then it is added first.
2161 * If q is in a column of the main tableau then the "+ q" can be accomplished
2162 * by setting the corresponding entry to the denominator of the constraint.
2163 * If q happens to be in a row of the main tableau, then the corresponding
2164 * row needs to be added instead (taking care of the denominators).
2165 * Note that this is very unlikely, but perhaps not entirely impossible.
2167 * The current value of the cut is known to be negative (or at least
2168 * non-positive), so row_sign is set accordingly.
2170 * Return the row of the cut or -1.
2172 static int add_parametric_cut(struct isl_tab *tab, int row,
2173 struct isl_context *context)
2175 struct isl_vec *div;
2176 int d;
2177 int i;
2178 int r;
2179 isl_int *r_row;
2180 int col;
2181 int n;
2182 unsigned off = 2 + tab->M;
2184 if (!context)
2185 return -1;
2187 div = get_row_parameter_div(tab, row);
2188 if (!div)
2189 return -1;
2191 n = tab->n_div - context->n_unknown;
2192 d = context->op->get_div(context, tab, div);
2193 isl_vec_free(div);
2194 if (d < 0)
2195 return -1;
2197 if (isl_tab_extend_cons(tab, 1) < 0)
2198 return -1;
2199 r = isl_tab_allocate_con(tab);
2200 if (r < 0)
2201 return -1;
2203 r_row = tab->mat->row[tab->con[r].index];
2204 isl_int_set(r_row[0], tab->mat->row[row][0]);
2205 isl_int_neg(r_row[1], tab->mat->row[row][1]);
2206 isl_int_fdiv_r(r_row[1], r_row[1], tab->mat->row[row][0]);
2207 isl_int_neg(r_row[1], r_row[1]);
2208 if (tab->M)
2209 isl_int_set_si(r_row[2], 0);
2210 for (i = 0; i < tab->n_param; ++i) {
2211 if (tab->var[i].is_row)
2212 continue;
2213 col = tab->var[i].index;
2214 isl_int_neg(r_row[off + col], tab->mat->row[row][off + col]);
2215 isl_int_fdiv_r(r_row[off + col], r_row[off + col],
2216 tab->mat->row[row][0]);
2217 isl_int_neg(r_row[off + col], r_row[off + col]);
2219 for (i = 0; i < tab->n_div; ++i) {
2220 if (tab->var[tab->n_var - tab->n_div + i].is_row)
2221 continue;
2222 col = tab->var[tab->n_var - tab->n_div + i].index;
2223 isl_int_neg(r_row[off + col], tab->mat->row[row][off + col]);
2224 isl_int_fdiv_r(r_row[off + col], r_row[off + col],
2225 tab->mat->row[row][0]);
2226 isl_int_neg(r_row[off + col], r_row[off + col]);
2228 for (i = 0; i < tab->n_col; ++i) {
2229 if (tab->col_var[i] >= 0 &&
2230 (tab->col_var[i] < tab->n_param ||
2231 tab->col_var[i] >= tab->n_var - tab->n_div))
2232 continue;
2233 isl_int_fdiv_r(r_row[off + i],
2234 tab->mat->row[row][off + i], tab->mat->row[row][0]);
2236 if (tab->var[tab->n_var - tab->n_div + d].is_row) {
2237 isl_int gcd;
2238 int d_row = tab->var[tab->n_var - tab->n_div + d].index;
2239 isl_int_init(gcd);
2240 isl_int_gcd(gcd, tab->mat->row[d_row][0], r_row[0]);
2241 isl_int_divexact(r_row[0], r_row[0], gcd);
2242 isl_int_divexact(gcd, tab->mat->row[d_row][0], gcd);
2243 isl_seq_combine(r_row + 1, gcd, r_row + 1,
2244 r_row[0], tab->mat->row[d_row] + 1,
2245 off - 1 + tab->n_col);
2246 isl_int_mul(r_row[0], r_row[0], tab->mat->row[d_row][0]);
2247 isl_int_clear(gcd);
2248 } else {
2249 col = tab->var[tab->n_var - tab->n_div + d].index;
2250 isl_int_set(r_row[off + col], tab->mat->row[row][0]);
2253 tab->con[r].is_nonneg = 1;
2254 if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r]) < 0)
2255 return -1;
2256 if (tab->row_sign)
2257 tab->row_sign[tab->con[r].index] = isl_tab_row_neg;
2259 row = tab->con[r].index;
2261 if (d >= n && context->op->detect_equalities(context, tab) < 0)
2262 return -1;
2264 return row;
2267 /* Construct a tableau for bmap that can be used for computing
2268 * the lexicographic minimum (or maximum) of bmap.
2269 * If not NULL, then dom is the domain where the minimum
2270 * should be computed. In this case, we set up a parametric
2271 * tableau with row signs (initialized to "unknown").
2272 * If M is set, then the tableau will use a big parameter.
2273 * If max is set, then a maximum should be computed instead of a minimum.
2274 * This means that for each variable x, the tableau will contain the variable
2275 * x' = M - x, rather than x' = M + x. This in turn means that the coefficient
2276 * of the variables in all constraints are negated prior to adding them
2277 * to the tableau.
2279 static __isl_give struct isl_tab *tab_for_lexmin(__isl_keep isl_basic_map *bmap,
2280 __isl_keep isl_basic_set *dom, unsigned M, int max)
2282 int i;
2283 struct isl_tab *tab;
2284 unsigned n_var;
2285 unsigned o_var;
2287 tab = isl_tab_alloc(bmap->ctx, 2 * bmap->n_eq + bmap->n_ineq + 1,
2288 isl_basic_map_total_dim(bmap), M);
2289 if (!tab)
2290 return NULL;
2292 tab->rational = ISL_F_ISSET(bmap, ISL_BASIC_MAP_RATIONAL);
2293 if (dom) {
2294 tab->n_param = isl_basic_set_total_dim(dom) - dom->n_div;
2295 tab->n_div = dom->n_div;
2296 tab->row_sign = isl_calloc_array(bmap->ctx,
2297 enum isl_tab_row_sign, tab->mat->n_row);
2298 if (tab->mat->n_row && !tab->row_sign)
2299 goto error;
2301 if (ISL_F_ISSET(bmap, ISL_BASIC_MAP_EMPTY)) {
2302 if (isl_tab_mark_empty(tab) < 0)
2303 goto error;
2304 return tab;
2307 for (i = tab->n_param; i < tab->n_var - tab->n_div; ++i) {
2308 tab->var[i].is_nonneg = 1;
2309 tab->var[i].frozen = 1;
2311 o_var = 1 + tab->n_param;
2312 n_var = tab->n_var - tab->n_param - tab->n_div;
2313 for (i = 0; i < bmap->n_eq; ++i) {
2314 if (max)
2315 isl_seq_neg(bmap->eq[i] + o_var,
2316 bmap->eq[i] + o_var, n_var);
2317 tab = add_lexmin_valid_eq(tab, bmap->eq[i]);
2318 if (max)
2319 isl_seq_neg(bmap->eq[i] + o_var,
2320 bmap->eq[i] + o_var, n_var);
2321 if (!tab || tab->empty)
2322 return tab;
2324 if (bmap->n_eq && restore_lexmin(tab) < 0)
2325 goto error;
2326 for (i = 0; i < bmap->n_ineq; ++i) {
2327 if (max)
2328 isl_seq_neg(bmap->ineq[i] + o_var,
2329 bmap->ineq[i] + o_var, n_var);
2330 tab = add_lexmin_ineq(tab, bmap->ineq[i]);
2331 if (max)
2332 isl_seq_neg(bmap->ineq[i] + o_var,
2333 bmap->ineq[i] + o_var, n_var);
2334 if (!tab || tab->empty)
2335 return tab;
2337 return tab;
2338 error:
2339 isl_tab_free(tab);
2340 return NULL;
2343 /* Given a main tableau where more than one row requires a split,
2344 * determine and return the "best" row to split on.
2346 * If any of the rows requiring a split only involves
2347 * variables that also appear in the context tableau,
2348 * then the negative part is guaranteed not to have a solution.
2349 * It is therefore best to split on any of these rows first.
2351 * Otherwise,
2352 * given two rows in the main tableau, if the inequality corresponding
2353 * to the first row is redundant with respect to that of the second row
2354 * in the current tableau, then it is better to split on the second row,
2355 * since in the positive part, both rows will be positive.
2356 * (In the negative part a pivot will have to be performed and just about
2357 * anything can happen to the sign of the other row.)
2359 * As a simple heuristic, we therefore select the row that makes the most
2360 * of the other rows redundant.
2362 * Perhaps it would also be useful to look at the number of constraints
2363 * that conflict with any given constraint.
2365 * best is the best row so far (-1 when we have not found any row yet).
2366 * best_r is the number of other rows made redundant by row best.
2367 * When best is still -1, bset_r is meaningless, but it is initialized
2368 * to some arbitrary value (0) anyway. Without this redundant initialization
2369 * valgrind may warn about uninitialized memory accesses when isl
2370 * is compiled with some versions of gcc.
2372 static int best_split(struct isl_tab *tab, struct isl_tab *context_tab)
2374 struct isl_tab_undo *snap;
2375 int split;
2376 int row;
2377 int best = -1;
2378 int best_r = 0;
2380 if (isl_tab_extend_cons(context_tab, 2) < 0)
2381 return -1;
2383 snap = isl_tab_snap(context_tab);
2385 for (split = tab->n_redundant; split < tab->n_row; ++split) {
2386 struct isl_tab_undo *snap2;
2387 struct isl_vec *ineq = NULL;
2388 int r = 0;
2389 int ok;
2391 if (!isl_tab_var_from_row(tab, split)->is_nonneg)
2392 continue;
2393 if (tab->row_sign[split] != isl_tab_row_any)
2394 continue;
2396 if (is_parametric_constant(tab, split))
2397 return split;
2399 ineq = get_row_parameter_ineq(tab, split);
2400 if (!ineq)
2401 return -1;
2402 ok = isl_tab_add_ineq(context_tab, ineq->el) >= 0;
2403 isl_vec_free(ineq);
2404 if (!ok)
2405 return -1;
2407 snap2 = isl_tab_snap(context_tab);
2409 for (row = tab->n_redundant; row < tab->n_row; ++row) {
2410 struct isl_tab_var *var;
2412 if (row == split)
2413 continue;
2414 if (!isl_tab_var_from_row(tab, row)->is_nonneg)
2415 continue;
2416 if (tab->row_sign[row] != isl_tab_row_any)
2417 continue;
2419 ineq = get_row_parameter_ineq(tab, row);
2420 if (!ineq)
2421 return -1;
2422 ok = isl_tab_add_ineq(context_tab, ineq->el) >= 0;
2423 isl_vec_free(ineq);
2424 if (!ok)
2425 return -1;
2426 var = &context_tab->con[context_tab->n_con - 1];
2427 if (!context_tab->empty &&
2428 !isl_tab_min_at_most_neg_one(context_tab, var))
2429 r++;
2430 if (isl_tab_rollback(context_tab, snap2) < 0)
2431 return -1;
2433 if (best == -1 || r > best_r) {
2434 best = split;
2435 best_r = r;
2437 if (isl_tab_rollback(context_tab, snap) < 0)
2438 return -1;
2441 return best;
2444 static struct isl_basic_set *context_lex_peek_basic_set(
2445 struct isl_context *context)
2447 struct isl_context_lex *clex = (struct isl_context_lex *)context;
2448 if (!clex->tab)
2449 return NULL;
2450 return isl_tab_peek_bset(clex->tab);
2453 static struct isl_tab *context_lex_peek_tab(struct isl_context *context)
2455 struct isl_context_lex *clex = (struct isl_context_lex *)context;
2456 return clex->tab;
2459 static void context_lex_add_eq(struct isl_context *context, isl_int *eq,
2460 int check, int update)
2462 struct isl_context_lex *clex = (struct isl_context_lex *)context;
2463 if (isl_tab_extend_cons(clex->tab, 2) < 0)
2464 goto error;
2465 if (add_lexmin_eq(clex->tab, eq) < 0)
2466 goto error;
2467 if (check) {
2468 int v = tab_has_valid_sample(clex->tab, eq, 1);
2469 if (v < 0)
2470 goto error;
2471 if (!v)
2472 clex->tab = check_integer_feasible(clex->tab);
2474 if (update)
2475 clex->tab = check_samples(clex->tab, eq, 1);
2476 return;
2477 error:
2478 isl_tab_free(clex->tab);
2479 clex->tab = NULL;
2482 static void context_lex_add_ineq(struct isl_context *context, isl_int *ineq,
2483 int check, int update)
2485 struct isl_context_lex *clex = (struct isl_context_lex *)context;
2486 if (isl_tab_extend_cons(clex->tab, 1) < 0)
2487 goto error;
2488 clex->tab = add_lexmin_ineq(clex->tab, ineq);
2489 if (check) {
2490 int v = tab_has_valid_sample(clex->tab, ineq, 0);
2491 if (v < 0)
2492 goto error;
2493 if (!v)
2494 clex->tab = check_integer_feasible(clex->tab);
2496 if (update)
2497 clex->tab = check_samples(clex->tab, ineq, 0);
2498 return;
2499 error:
2500 isl_tab_free(clex->tab);
2501 clex->tab = NULL;
2504 static isl_stat context_lex_add_ineq_wrap(void *user, isl_int *ineq)
2506 struct isl_context *context = (struct isl_context *)user;
2507 context_lex_add_ineq(context, ineq, 0, 0);
2508 return context->op->is_ok(context) ? isl_stat_ok : isl_stat_error;
2511 /* Check which signs can be obtained by "ineq" on all the currently
2512 * active sample values. See row_sign for more information.
2514 static enum isl_tab_row_sign tab_ineq_sign(struct isl_tab *tab, isl_int *ineq,
2515 int strict)
2517 int i;
2518 int sgn;
2519 isl_int tmp;
2520 enum isl_tab_row_sign res = isl_tab_row_unknown;
2522 isl_assert(tab->mat->ctx, tab->samples, return isl_tab_row_unknown);
2523 isl_assert(tab->mat->ctx, tab->samples->n_col == 1 + tab->n_var,
2524 return isl_tab_row_unknown);
2526 isl_int_init(tmp);
2527 for (i = tab->n_outside; i < tab->n_sample; ++i) {
2528 isl_seq_inner_product(tab->samples->row[i], ineq,
2529 1 + tab->n_var, &tmp);
2530 sgn = isl_int_sgn(tmp);
2531 if (sgn > 0 || (sgn == 0 && strict)) {
2532 if (res == isl_tab_row_unknown)
2533 res = isl_tab_row_pos;
2534 if (res == isl_tab_row_neg)
2535 res = isl_tab_row_any;
2537 if (sgn < 0) {
2538 if (res == isl_tab_row_unknown)
2539 res = isl_tab_row_neg;
2540 if (res == isl_tab_row_pos)
2541 res = isl_tab_row_any;
2543 if (res == isl_tab_row_any)
2544 break;
2546 isl_int_clear(tmp);
2548 return res;
2551 static enum isl_tab_row_sign context_lex_ineq_sign(struct isl_context *context,
2552 isl_int *ineq, int strict)
2554 struct isl_context_lex *clex = (struct isl_context_lex *)context;
2555 return tab_ineq_sign(clex->tab, ineq, strict);
2558 /* Check whether "ineq" can be added to the tableau without rendering
2559 * it infeasible.
2561 static int context_lex_test_ineq(struct isl_context *context, isl_int *ineq)
2563 struct isl_context_lex *clex = (struct isl_context_lex *)context;
2564 struct isl_tab_undo *snap;
2565 int feasible;
2567 if (!clex->tab)
2568 return -1;
2570 if (isl_tab_extend_cons(clex->tab, 1) < 0)
2571 return -1;
2573 snap = isl_tab_snap(clex->tab);
2574 if (isl_tab_push_basis(clex->tab) < 0)
2575 return -1;
2576 clex->tab = add_lexmin_ineq(clex->tab, ineq);
2577 clex->tab = check_integer_feasible(clex->tab);
2578 if (!clex->tab)
2579 return -1;
2580 feasible = !clex->tab->empty;
2581 if (isl_tab_rollback(clex->tab, snap) < 0)
2582 return -1;
2584 return feasible;
2587 static int context_lex_get_div(struct isl_context *context, struct isl_tab *tab,
2588 struct isl_vec *div)
2590 return get_div(tab, context, div);
2593 /* Insert a div specified by "div" to the context tableau at position "pos" and
2594 * return isl_bool_true if the div is obviously non-negative.
2595 * context_tab_add_div will always return isl_bool_true, because all variables
2596 * in a isl_context_lex tableau are non-negative.
2597 * However, if we are using a big parameter in the context, then this only
2598 * reflects the non-negativity of the variable used to _encode_ the
2599 * div, i.e., div' = M + div, so we can't draw any conclusions.
2601 static isl_bool context_lex_insert_div(struct isl_context *context, int pos,
2602 __isl_keep isl_vec *div)
2604 struct isl_context_lex *clex = (struct isl_context_lex *)context;
2605 isl_bool nonneg;
2606 nonneg = context_tab_insert_div(clex->tab, pos, div,
2607 context_lex_add_ineq_wrap, context);
2608 if (nonneg < 0)
2609 return isl_bool_error;
2610 if (clex->tab->M)
2611 return isl_bool_false;
2612 return nonneg;
2615 static int context_lex_detect_equalities(struct isl_context *context,
2616 struct isl_tab *tab)
2618 return 0;
2621 static int context_lex_best_split(struct isl_context *context,
2622 struct isl_tab *tab)
2624 struct isl_context_lex *clex = (struct isl_context_lex *)context;
2625 struct isl_tab_undo *snap;
2626 int r;
2628 snap = isl_tab_snap(clex->tab);
2629 if (isl_tab_push_basis(clex->tab) < 0)
2630 return -1;
2631 r = best_split(tab, clex->tab);
2633 if (r >= 0 && isl_tab_rollback(clex->tab, snap) < 0)
2634 return -1;
2636 return r;
2639 static int context_lex_is_empty(struct isl_context *context)
2641 struct isl_context_lex *clex = (struct isl_context_lex *)context;
2642 if (!clex->tab)
2643 return -1;
2644 return clex->tab->empty;
2647 static void *context_lex_save(struct isl_context *context)
2649 struct isl_context_lex *clex = (struct isl_context_lex *)context;
2650 struct isl_tab_undo *snap;
2652 snap = isl_tab_snap(clex->tab);
2653 if (isl_tab_push_basis(clex->tab) < 0)
2654 return NULL;
2655 if (isl_tab_save_samples(clex->tab) < 0)
2656 return NULL;
2658 return snap;
2661 static void context_lex_restore(struct isl_context *context, void *save)
2663 struct isl_context_lex *clex = (struct isl_context_lex *)context;
2664 if (isl_tab_rollback(clex->tab, (struct isl_tab_undo *)save) < 0) {
2665 isl_tab_free(clex->tab);
2666 clex->tab = NULL;
2670 static void context_lex_discard(void *save)
2674 static int context_lex_is_ok(struct isl_context *context)
2676 struct isl_context_lex *clex = (struct isl_context_lex *)context;
2677 return !!clex->tab;
2680 /* For each variable in the context tableau, check if the variable can
2681 * only attain non-negative values. If so, mark the parameter as non-negative
2682 * in the main tableau. This allows for a more direct identification of some
2683 * cases of violated constraints.
2685 static struct isl_tab *tab_detect_nonnegative_parameters(struct isl_tab *tab,
2686 struct isl_tab *context_tab)
2688 int i;
2689 struct isl_tab_undo *snap;
2690 struct isl_vec *ineq = NULL;
2691 struct isl_tab_var *var;
2692 int n;
2694 if (context_tab->n_var == 0)
2695 return tab;
2697 ineq = isl_vec_alloc(tab->mat->ctx, 1 + context_tab->n_var);
2698 if (!ineq)
2699 goto error;
2701 if (isl_tab_extend_cons(context_tab, 1) < 0)
2702 goto error;
2704 snap = isl_tab_snap(context_tab);
2706 n = 0;
2707 isl_seq_clr(ineq->el, ineq->size);
2708 for (i = 0; i < context_tab->n_var; ++i) {
2709 isl_int_set_si(ineq->el[1 + i], 1);
2710 if (isl_tab_add_ineq(context_tab, ineq->el) < 0)
2711 goto error;
2712 var = &context_tab->con[context_tab->n_con - 1];
2713 if (!context_tab->empty &&
2714 !isl_tab_min_at_most_neg_one(context_tab, var)) {
2715 int j = i;
2716 if (i >= tab->n_param)
2717 j = i - tab->n_param + tab->n_var - tab->n_div;
2718 tab->var[j].is_nonneg = 1;
2719 n++;
2721 isl_int_set_si(ineq->el[1 + i], 0);
2722 if (isl_tab_rollback(context_tab, snap) < 0)
2723 goto error;
2726 if (context_tab->M && n == context_tab->n_var) {
2727 context_tab->mat = isl_mat_drop_cols(context_tab->mat, 2, 1);
2728 context_tab->M = 0;
2731 isl_vec_free(ineq);
2732 return tab;
2733 error:
2734 isl_vec_free(ineq);
2735 isl_tab_free(tab);
2736 return NULL;
2739 static struct isl_tab *context_lex_detect_nonnegative_parameters(
2740 struct isl_context *context, struct isl_tab *tab)
2742 struct isl_context_lex *clex = (struct isl_context_lex *)context;
2743 struct isl_tab_undo *snap;
2745 if (!tab)
2746 return NULL;
2748 snap = isl_tab_snap(clex->tab);
2749 if (isl_tab_push_basis(clex->tab) < 0)
2750 goto error;
2752 tab = tab_detect_nonnegative_parameters(tab, clex->tab);
2754 if (isl_tab_rollback(clex->tab, snap) < 0)
2755 goto error;
2757 return tab;
2758 error:
2759 isl_tab_free(tab);
2760 return NULL;
2763 static void context_lex_invalidate(struct isl_context *context)
2765 struct isl_context_lex *clex = (struct isl_context_lex *)context;
2766 isl_tab_free(clex->tab);
2767 clex->tab = NULL;
2770 static __isl_null struct isl_context *context_lex_free(
2771 struct isl_context *context)
2773 struct isl_context_lex *clex = (struct isl_context_lex *)context;
2774 isl_tab_free(clex->tab);
2775 free(clex);
2777 return NULL;
2780 struct isl_context_op isl_context_lex_op = {
2781 context_lex_detect_nonnegative_parameters,
2782 context_lex_peek_basic_set,
2783 context_lex_peek_tab,
2784 context_lex_add_eq,
2785 context_lex_add_ineq,
2786 context_lex_ineq_sign,
2787 context_lex_test_ineq,
2788 context_lex_get_div,
2789 context_lex_insert_div,
2790 context_lex_detect_equalities,
2791 context_lex_best_split,
2792 context_lex_is_empty,
2793 context_lex_is_ok,
2794 context_lex_save,
2795 context_lex_restore,
2796 context_lex_discard,
2797 context_lex_invalidate,
2798 context_lex_free,
2801 static struct isl_tab *context_tab_for_lexmin(__isl_take isl_basic_set *bset)
2803 struct isl_tab *tab;
2805 if (!bset)
2806 return NULL;
2807 tab = tab_for_lexmin(bset_to_bmap(bset), NULL, 1, 0);
2808 if (isl_tab_track_bset(tab, bset) < 0)
2809 goto error;
2810 tab = isl_tab_init_samples(tab);
2811 return tab;
2812 error:
2813 isl_tab_free(tab);
2814 return NULL;
2817 static struct isl_context *isl_context_lex_alloc(struct isl_basic_set *dom)
2819 struct isl_context_lex *clex;
2821 if (!dom)
2822 return NULL;
2824 clex = isl_alloc_type(dom->ctx, struct isl_context_lex);
2825 if (!clex)
2826 return NULL;
2828 clex->context.op = &isl_context_lex_op;
2830 clex->tab = context_tab_for_lexmin(isl_basic_set_copy(dom));
2831 if (restore_lexmin(clex->tab) < 0)
2832 goto error;
2833 clex->tab = check_integer_feasible(clex->tab);
2834 if (!clex->tab)
2835 goto error;
2837 return &clex->context;
2838 error:
2839 clex->context.op->free(&clex->context);
2840 return NULL;
2843 /* Representation of the context when using generalized basis reduction.
2845 * "shifted" contains the offsets of the unit hypercubes that lie inside the
2846 * context. Any rational point in "shifted" can therefore be rounded
2847 * up to an integer point in the context.
2848 * If the context is constrained by any equality, then "shifted" is not used
2849 * as it would be empty.
2851 struct isl_context_gbr {
2852 struct isl_context context;
2853 struct isl_tab *tab;
2854 struct isl_tab *shifted;
2855 struct isl_tab *cone;
2858 static struct isl_tab *context_gbr_detect_nonnegative_parameters(
2859 struct isl_context *context, struct isl_tab *tab)
2861 struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
2862 if (!tab)
2863 return NULL;
2864 return tab_detect_nonnegative_parameters(tab, cgbr->tab);
2867 static struct isl_basic_set *context_gbr_peek_basic_set(
2868 struct isl_context *context)
2870 struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
2871 if (!cgbr->tab)
2872 return NULL;
2873 return isl_tab_peek_bset(cgbr->tab);
2876 static struct isl_tab *context_gbr_peek_tab(struct isl_context *context)
2878 struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
2879 return cgbr->tab;
2882 /* Initialize the "shifted" tableau of the context, which
2883 * contains the constraints of the original tableau shifted
2884 * by the sum of all negative coefficients. This ensures
2885 * that any rational point in the shifted tableau can
2886 * be rounded up to yield an integer point in the original tableau.
2888 static void gbr_init_shifted(struct isl_context_gbr *cgbr)
2890 int i, j;
2891 struct isl_vec *cst;
2892 struct isl_basic_set *bset = isl_tab_peek_bset(cgbr->tab);
2893 unsigned dim = isl_basic_set_total_dim(bset);
2895 cst = isl_vec_alloc(cgbr->tab->mat->ctx, bset->n_ineq);
2896 if (!cst)
2897 return;
2899 for (i = 0; i < bset->n_ineq; ++i) {
2900 isl_int_set(cst->el[i], bset->ineq[i][0]);
2901 for (j = 0; j < dim; ++j) {
2902 if (!isl_int_is_neg(bset->ineq[i][1 + j]))
2903 continue;
2904 isl_int_add(bset->ineq[i][0], bset->ineq[i][0],
2905 bset->ineq[i][1 + j]);
2909 cgbr->shifted = isl_tab_from_basic_set(bset, 0);
2911 for (i = 0; i < bset->n_ineq; ++i)
2912 isl_int_set(bset->ineq[i][0], cst->el[i]);
2914 isl_vec_free(cst);
2917 /* Check if the shifted tableau is non-empty, and if so
2918 * use the sample point to construct an integer point
2919 * of the context tableau.
2921 static struct isl_vec *gbr_get_shifted_sample(struct isl_context_gbr *cgbr)
2923 struct isl_vec *sample;
2925 if (!cgbr->shifted)
2926 gbr_init_shifted(cgbr);
2927 if (!cgbr->shifted)
2928 return NULL;
2929 if (cgbr->shifted->empty)
2930 return isl_vec_alloc(cgbr->tab->mat->ctx, 0);
2932 sample = isl_tab_get_sample_value(cgbr->shifted);
2933 sample = isl_vec_ceil(sample);
2935 return sample;
2938 static __isl_give isl_basic_set *drop_constant_terms(
2939 __isl_take isl_basic_set *bset)
2941 int i;
2943 if (!bset)
2944 return NULL;
2946 for (i = 0; i < bset->n_eq; ++i)
2947 isl_int_set_si(bset->eq[i][0], 0);
2949 for (i = 0; i < bset->n_ineq; ++i)
2950 isl_int_set_si(bset->ineq[i][0], 0);
2952 return bset;
2955 static int use_shifted(struct isl_context_gbr *cgbr)
2957 if (!cgbr->tab)
2958 return 0;
2959 return cgbr->tab->bmap->n_eq == 0 && cgbr->tab->bmap->n_div == 0;
2962 static struct isl_vec *gbr_get_sample(struct isl_context_gbr *cgbr)
2964 struct isl_basic_set *bset;
2965 struct isl_basic_set *cone;
2967 if (isl_tab_sample_is_integer(cgbr->tab))
2968 return isl_tab_get_sample_value(cgbr->tab);
2970 if (use_shifted(cgbr)) {
2971 struct isl_vec *sample;
2973 sample = gbr_get_shifted_sample(cgbr);
2974 if (!sample || sample->size > 0)
2975 return sample;
2977 isl_vec_free(sample);
2980 if (!cgbr->cone) {
2981 bset = isl_tab_peek_bset(cgbr->tab);
2982 cgbr->cone = isl_tab_from_recession_cone(bset, 0);
2983 if (!cgbr->cone)
2984 return NULL;
2985 if (isl_tab_track_bset(cgbr->cone,
2986 isl_basic_set_copy(bset)) < 0)
2987 return NULL;
2989 if (isl_tab_detect_implicit_equalities(cgbr->cone) < 0)
2990 return NULL;
2992 if (cgbr->cone->n_dead == cgbr->cone->n_col) {
2993 struct isl_vec *sample;
2994 struct isl_tab_undo *snap;
2996 if (cgbr->tab->basis) {
2997 if (cgbr->tab->basis->n_col != 1 + cgbr->tab->n_var) {
2998 isl_mat_free(cgbr->tab->basis);
2999 cgbr->tab->basis = NULL;
3001 cgbr->tab->n_zero = 0;
3002 cgbr->tab->n_unbounded = 0;
3005 snap = isl_tab_snap(cgbr->tab);
3007 sample = isl_tab_sample(cgbr->tab);
3009 if (!sample || isl_tab_rollback(cgbr->tab, snap) < 0) {
3010 isl_vec_free(sample);
3011 return NULL;
3014 return sample;
3017 cone = isl_basic_set_dup(isl_tab_peek_bset(cgbr->cone));
3018 cone = drop_constant_terms(cone);
3019 cone = isl_basic_set_update_from_tab(cone, cgbr->cone);
3020 cone = isl_basic_set_underlying_set(cone);
3021 cone = isl_basic_set_gauss(cone, NULL);
3023 bset = isl_basic_set_dup(isl_tab_peek_bset(cgbr->tab));
3024 bset = isl_basic_set_update_from_tab(bset, cgbr->tab);
3025 bset = isl_basic_set_underlying_set(bset);
3026 bset = isl_basic_set_gauss(bset, NULL);
3028 return isl_basic_set_sample_with_cone(bset, cone);
3031 static void check_gbr_integer_feasible(struct isl_context_gbr *cgbr)
3033 struct isl_vec *sample;
3035 if (!cgbr->tab)
3036 return;
3038 if (cgbr->tab->empty)
3039 return;
3041 sample = gbr_get_sample(cgbr);
3042 if (!sample)
3043 goto error;
3045 if (sample->size == 0) {
3046 isl_vec_free(sample);
3047 if (isl_tab_mark_empty(cgbr->tab) < 0)
3048 goto error;
3049 return;
3052 if (isl_tab_add_sample(cgbr->tab, sample) < 0)
3053 goto error;
3055 return;
3056 error:
3057 isl_tab_free(cgbr->tab);
3058 cgbr->tab = NULL;
3061 static struct isl_tab *add_gbr_eq(struct isl_tab *tab, isl_int *eq)
3063 if (!tab)
3064 return NULL;
3066 if (isl_tab_extend_cons(tab, 2) < 0)
3067 goto error;
3069 if (isl_tab_add_eq(tab, eq) < 0)
3070 goto error;
3072 return tab;
3073 error:
3074 isl_tab_free(tab);
3075 return NULL;
3078 /* Add the equality described by "eq" to the context.
3079 * If "check" is set, then we check if the context is empty after
3080 * adding the equality.
3081 * If "update" is set, then we check if the samples are still valid.
3083 * We do not explicitly add shifted copies of the equality to
3084 * cgbr->shifted since they would conflict with each other.
3085 * Instead, we directly mark cgbr->shifted empty.
3087 static void context_gbr_add_eq(struct isl_context *context, isl_int *eq,
3088 int check, int update)
3090 struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3092 cgbr->tab = add_gbr_eq(cgbr->tab, eq);
3094 if (cgbr->shifted && !cgbr->shifted->empty && use_shifted(cgbr)) {
3095 if (isl_tab_mark_empty(cgbr->shifted) < 0)
3096 goto error;
3099 if (cgbr->cone && cgbr->cone->n_col != cgbr->cone->n_dead) {
3100 if (isl_tab_extend_cons(cgbr->cone, 2) < 0)
3101 goto error;
3102 if (isl_tab_add_eq(cgbr->cone, eq) < 0)
3103 goto error;
3106 if (check) {
3107 int v = tab_has_valid_sample(cgbr->tab, eq, 1);
3108 if (v < 0)
3109 goto error;
3110 if (!v)
3111 check_gbr_integer_feasible(cgbr);
3113 if (update)
3114 cgbr->tab = check_samples(cgbr->tab, eq, 1);
3115 return;
3116 error:
3117 isl_tab_free(cgbr->tab);
3118 cgbr->tab = NULL;
3121 static void add_gbr_ineq(struct isl_context_gbr *cgbr, isl_int *ineq)
3123 if (!cgbr->tab)
3124 return;
3126 if (isl_tab_extend_cons(cgbr->tab, 1) < 0)
3127 goto error;
3129 if (isl_tab_add_ineq(cgbr->tab, ineq) < 0)
3130 goto error;
3132 if (cgbr->shifted && !cgbr->shifted->empty && use_shifted(cgbr)) {
3133 int i;
3134 unsigned dim;
3135 dim = isl_basic_map_total_dim(cgbr->tab->bmap);
3137 if (isl_tab_extend_cons(cgbr->shifted, 1) < 0)
3138 goto error;
3140 for (i = 0; i < dim; ++i) {
3141 if (!isl_int_is_neg(ineq[1 + i]))
3142 continue;
3143 isl_int_add(ineq[0], ineq[0], ineq[1 + i]);
3146 if (isl_tab_add_ineq(cgbr->shifted, ineq) < 0)
3147 goto error;
3149 for (i = 0; i < dim; ++i) {
3150 if (!isl_int_is_neg(ineq[1 + i]))
3151 continue;
3152 isl_int_sub(ineq[0], ineq[0], ineq[1 + i]);
3156 if (cgbr->cone && cgbr->cone->n_col != cgbr->cone->n_dead) {
3157 if (isl_tab_extend_cons(cgbr->cone, 1) < 0)
3158 goto error;
3159 if (isl_tab_add_ineq(cgbr->cone, ineq) < 0)
3160 goto error;
3163 return;
3164 error:
3165 isl_tab_free(cgbr->tab);
3166 cgbr->tab = NULL;
3169 static void context_gbr_add_ineq(struct isl_context *context, isl_int *ineq,
3170 int check, int update)
3172 struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3174 add_gbr_ineq(cgbr, ineq);
3175 if (!cgbr->tab)
3176 return;
3178 if (check) {
3179 int v = tab_has_valid_sample(cgbr->tab, ineq, 0);
3180 if (v < 0)
3181 goto error;
3182 if (!v)
3183 check_gbr_integer_feasible(cgbr);
3185 if (update)
3186 cgbr->tab = check_samples(cgbr->tab, ineq, 0);
3187 return;
3188 error:
3189 isl_tab_free(cgbr->tab);
3190 cgbr->tab = NULL;
3193 static isl_stat context_gbr_add_ineq_wrap(void *user, isl_int *ineq)
3195 struct isl_context *context = (struct isl_context *)user;
3196 context_gbr_add_ineq(context, ineq, 0, 0);
3197 return context->op->is_ok(context) ? isl_stat_ok : isl_stat_error;
3200 static enum isl_tab_row_sign context_gbr_ineq_sign(struct isl_context *context,
3201 isl_int *ineq, int strict)
3203 struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3204 return tab_ineq_sign(cgbr->tab, ineq, strict);
3207 /* Check whether "ineq" can be added to the tableau without rendering
3208 * it infeasible.
3210 static int context_gbr_test_ineq(struct isl_context *context, isl_int *ineq)
3212 struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3213 struct isl_tab_undo *snap;
3214 struct isl_tab_undo *shifted_snap = NULL;
3215 struct isl_tab_undo *cone_snap = NULL;
3216 int feasible;
3218 if (!cgbr->tab)
3219 return -1;
3221 if (isl_tab_extend_cons(cgbr->tab, 1) < 0)
3222 return -1;
3224 snap = isl_tab_snap(cgbr->tab);
3225 if (cgbr->shifted)
3226 shifted_snap = isl_tab_snap(cgbr->shifted);
3227 if (cgbr->cone)
3228 cone_snap = isl_tab_snap(cgbr->cone);
3229 add_gbr_ineq(cgbr, ineq);
3230 check_gbr_integer_feasible(cgbr);
3231 if (!cgbr->tab)
3232 return -1;
3233 feasible = !cgbr->tab->empty;
3234 if (isl_tab_rollback(cgbr->tab, snap) < 0)
3235 return -1;
3236 if (shifted_snap) {
3237 if (isl_tab_rollback(cgbr->shifted, shifted_snap))
3238 return -1;
3239 } else if (cgbr->shifted) {
3240 isl_tab_free(cgbr->shifted);
3241 cgbr->shifted = NULL;
3243 if (cone_snap) {
3244 if (isl_tab_rollback(cgbr->cone, cone_snap))
3245 return -1;
3246 } else if (cgbr->cone) {
3247 isl_tab_free(cgbr->cone);
3248 cgbr->cone = NULL;
3251 return feasible;
3254 /* Return the column of the last of the variables associated to
3255 * a column that has a non-zero coefficient.
3256 * This function is called in a context where only coefficients
3257 * of parameters or divs can be non-zero.
3259 static int last_non_zero_var_col(struct isl_tab *tab, isl_int *p)
3261 int i;
3262 int col;
3264 if (tab->n_var == 0)
3265 return -1;
3267 for (i = tab->n_var - 1; i >= 0; --i) {
3268 if (i >= tab->n_param && i < tab->n_var - tab->n_div)
3269 continue;
3270 if (tab->var[i].is_row)
3271 continue;
3272 col = tab->var[i].index;
3273 if (!isl_int_is_zero(p[col]))
3274 return col;
3277 return -1;
3280 /* Look through all the recently added equalities in the context
3281 * to see if we can propagate any of them to the main tableau.
3283 * The newly added equalities in the context are encoded as pairs
3284 * of inequalities starting at inequality "first".
3286 * We tentatively add each of these equalities to the main tableau
3287 * and if this happens to result in a row with a final coefficient
3288 * that is one or negative one, we use it to kill a column
3289 * in the main tableau. Otherwise, we discard the tentatively
3290 * added row.
3291 * This tentative addition of equality constraints turns
3292 * on the undo facility of the tableau. Turn it off again
3293 * at the end, assuming it was turned off to begin with.
3295 * Return 0 on success and -1 on failure.
3297 static int propagate_equalities(struct isl_context_gbr *cgbr,
3298 struct isl_tab *tab, unsigned first)
3300 int i;
3301 struct isl_vec *eq = NULL;
3302 isl_bool needs_undo;
3304 needs_undo = isl_tab_need_undo(tab);
3305 if (needs_undo < 0)
3306 goto error;
3307 eq = isl_vec_alloc(tab->mat->ctx, 1 + tab->n_var);
3308 if (!eq)
3309 goto error;
3311 if (isl_tab_extend_cons(tab, (cgbr->tab->bmap->n_ineq - first)/2) < 0)
3312 goto error;
3314 isl_seq_clr(eq->el + 1 + tab->n_param,
3315 tab->n_var - tab->n_param - tab->n_div);
3316 for (i = first; i < cgbr->tab->bmap->n_ineq; i += 2) {
3317 int j;
3318 int r;
3319 struct isl_tab_undo *snap;
3320 snap = isl_tab_snap(tab);
3322 isl_seq_cpy(eq->el, cgbr->tab->bmap->ineq[i], 1 + tab->n_param);
3323 isl_seq_cpy(eq->el + 1 + tab->n_var - tab->n_div,
3324 cgbr->tab->bmap->ineq[i] + 1 + tab->n_param,
3325 tab->n_div);
3327 r = isl_tab_add_row(tab, eq->el);
3328 if (r < 0)
3329 goto error;
3330 r = tab->con[r].index;
3331 j = last_non_zero_var_col(tab, tab->mat->row[r] + 2 + tab->M);
3332 if (j < 0 || j < tab->n_dead ||
3333 !isl_int_is_one(tab->mat->row[r][0]) ||
3334 (!isl_int_is_one(tab->mat->row[r][2 + tab->M + j]) &&
3335 !isl_int_is_negone(tab->mat->row[r][2 + tab->M + j]))) {
3336 if (isl_tab_rollback(tab, snap) < 0)
3337 goto error;
3338 continue;
3340 if (isl_tab_pivot(tab, r, j) < 0)
3341 goto error;
3342 if (isl_tab_kill_col(tab, j) < 0)
3343 goto error;
3345 if (restore_lexmin(tab) < 0)
3346 goto error;
3349 if (!needs_undo)
3350 isl_tab_clear_undo(tab);
3351 isl_vec_free(eq);
3353 return 0;
3354 error:
3355 isl_vec_free(eq);
3356 isl_tab_free(cgbr->tab);
3357 cgbr->tab = NULL;
3358 return -1;
3361 static int context_gbr_detect_equalities(struct isl_context *context,
3362 struct isl_tab *tab)
3364 struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3365 unsigned n_ineq;
3367 if (!cgbr->cone) {
3368 struct isl_basic_set *bset = isl_tab_peek_bset(cgbr->tab);
3369 cgbr->cone = isl_tab_from_recession_cone(bset, 0);
3370 if (!cgbr->cone)
3371 goto error;
3372 if (isl_tab_track_bset(cgbr->cone,
3373 isl_basic_set_copy(bset)) < 0)
3374 goto error;
3376 if (isl_tab_detect_implicit_equalities(cgbr->cone) < 0)
3377 goto error;
3379 n_ineq = cgbr->tab->bmap->n_ineq;
3380 cgbr->tab = isl_tab_detect_equalities(cgbr->tab, cgbr->cone);
3381 if (!cgbr->tab)
3382 return -1;
3383 if (cgbr->tab->bmap->n_ineq > n_ineq &&
3384 propagate_equalities(cgbr, tab, n_ineq) < 0)
3385 return -1;
3387 return 0;
3388 error:
3389 isl_tab_free(cgbr->tab);
3390 cgbr->tab = NULL;
3391 return -1;
3394 static int context_gbr_get_div(struct isl_context *context, struct isl_tab *tab,
3395 struct isl_vec *div)
3397 return get_div(tab, context, div);
3400 static isl_bool context_gbr_insert_div(struct isl_context *context, int pos,
3401 __isl_keep isl_vec *div)
3403 struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3404 if (cgbr->cone) {
3405 int r, n_div, o_div;
3407 n_div = isl_basic_map_dim(cgbr->cone->bmap, isl_dim_div);
3408 o_div = cgbr->cone->n_var - n_div;
3410 if (isl_tab_extend_cons(cgbr->cone, 3) < 0)
3411 return isl_bool_error;
3412 if (isl_tab_extend_vars(cgbr->cone, 1) < 0)
3413 return isl_bool_error;
3414 if ((r = isl_tab_insert_var(cgbr->cone, pos)) <0)
3415 return isl_bool_error;
3417 cgbr->cone->bmap = isl_basic_map_insert_div(cgbr->cone->bmap,
3418 r - o_div, div);
3419 if (!cgbr->cone->bmap)
3420 return isl_bool_error;
3421 if (isl_tab_push_var(cgbr->cone, isl_tab_undo_bmap_div,
3422 &cgbr->cone->var[r]) < 0)
3423 return isl_bool_error;
3425 return context_tab_insert_div(cgbr->tab, pos, div,
3426 context_gbr_add_ineq_wrap, context);
3429 static int context_gbr_best_split(struct isl_context *context,
3430 struct isl_tab *tab)
3432 struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3433 struct isl_tab_undo *snap;
3434 int r;
3436 snap = isl_tab_snap(cgbr->tab);
3437 r = best_split(tab, cgbr->tab);
3439 if (r >= 0 && isl_tab_rollback(cgbr->tab, snap) < 0)
3440 return -1;
3442 return r;
3445 static int context_gbr_is_empty(struct isl_context *context)
3447 struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3448 if (!cgbr->tab)
3449 return -1;
3450 return cgbr->tab->empty;
3453 struct isl_gbr_tab_undo {
3454 struct isl_tab_undo *tab_snap;
3455 struct isl_tab_undo *shifted_snap;
3456 struct isl_tab_undo *cone_snap;
3459 static void *context_gbr_save(struct isl_context *context)
3461 struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3462 struct isl_gbr_tab_undo *snap;
3464 if (!cgbr->tab)
3465 return NULL;
3467 snap = isl_alloc_type(cgbr->tab->mat->ctx, struct isl_gbr_tab_undo);
3468 if (!snap)
3469 return NULL;
3471 snap->tab_snap = isl_tab_snap(cgbr->tab);
3472 if (isl_tab_save_samples(cgbr->tab) < 0)
3473 goto error;
3475 if (cgbr->shifted)
3476 snap->shifted_snap = isl_tab_snap(cgbr->shifted);
3477 else
3478 snap->shifted_snap = NULL;
3480 if (cgbr->cone)
3481 snap->cone_snap = isl_tab_snap(cgbr->cone);
3482 else
3483 snap->cone_snap = NULL;
3485 return snap;
3486 error:
3487 free(snap);
3488 return NULL;
3491 static void context_gbr_restore(struct isl_context *context, void *save)
3493 struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3494 struct isl_gbr_tab_undo *snap = (struct isl_gbr_tab_undo *)save;
3495 if (!snap)
3496 goto error;
3497 if (isl_tab_rollback(cgbr->tab, snap->tab_snap) < 0)
3498 goto error;
3500 if (snap->shifted_snap) {
3501 if (isl_tab_rollback(cgbr->shifted, snap->shifted_snap) < 0)
3502 goto error;
3503 } else if (cgbr->shifted) {
3504 isl_tab_free(cgbr->shifted);
3505 cgbr->shifted = NULL;
3508 if (snap->cone_snap) {
3509 if (isl_tab_rollback(cgbr->cone, snap->cone_snap) < 0)
3510 goto error;
3511 } else if (cgbr->cone) {
3512 isl_tab_free(cgbr->cone);
3513 cgbr->cone = NULL;
3516 free(snap);
3518 return;
3519 error:
3520 free(snap);
3521 isl_tab_free(cgbr->tab);
3522 cgbr->tab = NULL;
3525 static void context_gbr_discard(void *save)
3527 struct isl_gbr_tab_undo *snap = (struct isl_gbr_tab_undo *)save;
3528 free(snap);
3531 static int context_gbr_is_ok(struct isl_context *context)
3533 struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3534 return !!cgbr->tab;
3537 static void context_gbr_invalidate(struct isl_context *context)
3539 struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3540 isl_tab_free(cgbr->tab);
3541 cgbr->tab = NULL;
3544 static __isl_null struct isl_context *context_gbr_free(
3545 struct isl_context *context)
3547 struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3548 isl_tab_free(cgbr->tab);
3549 isl_tab_free(cgbr->shifted);
3550 isl_tab_free(cgbr->cone);
3551 free(cgbr);
3553 return NULL;
3556 struct isl_context_op isl_context_gbr_op = {
3557 context_gbr_detect_nonnegative_parameters,
3558 context_gbr_peek_basic_set,
3559 context_gbr_peek_tab,
3560 context_gbr_add_eq,
3561 context_gbr_add_ineq,
3562 context_gbr_ineq_sign,
3563 context_gbr_test_ineq,
3564 context_gbr_get_div,
3565 context_gbr_insert_div,
3566 context_gbr_detect_equalities,
3567 context_gbr_best_split,
3568 context_gbr_is_empty,
3569 context_gbr_is_ok,
3570 context_gbr_save,
3571 context_gbr_restore,
3572 context_gbr_discard,
3573 context_gbr_invalidate,
3574 context_gbr_free,
3577 static struct isl_context *isl_context_gbr_alloc(__isl_keep isl_basic_set *dom)
3579 struct isl_context_gbr *cgbr;
3581 if (!dom)
3582 return NULL;
3584 cgbr = isl_calloc_type(dom->ctx, struct isl_context_gbr);
3585 if (!cgbr)
3586 return NULL;
3588 cgbr->context.op = &isl_context_gbr_op;
3590 cgbr->shifted = NULL;
3591 cgbr->cone = NULL;
3592 cgbr->tab = isl_tab_from_basic_set(dom, 1);
3593 cgbr->tab = isl_tab_init_samples(cgbr->tab);
3594 if (!cgbr->tab)
3595 goto error;
3596 check_gbr_integer_feasible(cgbr);
3598 return &cgbr->context;
3599 error:
3600 cgbr->context.op->free(&cgbr->context);
3601 return NULL;
3604 /* Allocate a context corresponding to "dom".
3605 * The representation specific fields are initialized by
3606 * isl_context_lex_alloc or isl_context_gbr_alloc.
3607 * The shared "n_unknown" field is initialized to the number
3608 * of final unknown integer divisions in "dom".
3610 static struct isl_context *isl_context_alloc(__isl_keep isl_basic_set *dom)
3612 struct isl_context *context;
3613 int first;
3615 if (!dom)
3616 return NULL;
3618 if (dom->ctx->opt->context == ISL_CONTEXT_LEXMIN)
3619 context = isl_context_lex_alloc(dom);
3620 else
3621 context = isl_context_gbr_alloc(dom);
3623 if (!context)
3624 return NULL;
3626 first = isl_basic_set_first_unknown_div(dom);
3627 if (first < 0)
3628 return context->op->free(context);
3629 context->n_unknown = isl_basic_set_dim(dom, isl_dim_div) - first;
3631 return context;
3634 /* Initialize some common fields of "sol", which keeps track
3635 * of the solution of an optimization problem on "bmap" over
3636 * the domain "dom".
3637 * If "max" is set, then a maximization problem is being solved, rather than
3638 * a minimization problem, which means that the variables in the
3639 * tableau have value "M - x" rather than "M + x".
3641 static isl_stat sol_init(struct isl_sol *sol, __isl_keep isl_basic_map *bmap,
3642 __isl_keep isl_basic_set *dom, int max)
3644 sol->rational = ISL_F_ISSET(bmap, ISL_BASIC_MAP_RATIONAL);
3645 sol->dec_level.callback.run = &sol_dec_level_wrap;
3646 sol->dec_level.sol = sol;
3647 sol->max = max;
3648 sol->n_out = isl_basic_map_dim(bmap, isl_dim_out);
3649 sol->space = isl_basic_map_get_space(bmap);
3651 sol->context = isl_context_alloc(dom);
3652 if (!sol->space || !sol->context)
3653 return isl_stat_error;
3655 return isl_stat_ok;
3658 /* Construct an isl_sol_map structure for accumulating the solution.
3659 * If track_empty is set, then we also keep track of the parts
3660 * of the context where there is no solution.
3661 * If max is set, then we are solving a maximization, rather than
3662 * a minimization problem, which means that the variables in the
3663 * tableau have value "M - x" rather than "M + x".
3665 static struct isl_sol *sol_map_init(__isl_keep isl_basic_map *bmap,
3666 __isl_take isl_basic_set *dom, int track_empty, int max)
3668 struct isl_sol_map *sol_map = NULL;
3669 isl_space *space;
3671 if (!bmap)
3672 goto error;
3674 sol_map = isl_calloc_type(bmap->ctx, struct isl_sol_map);
3675 if (!sol_map)
3676 goto error;
3678 sol_map->sol.free = &sol_map_free;
3679 if (sol_init(&sol_map->sol, bmap, dom, max) < 0)
3680 goto error;
3681 sol_map->sol.add = &sol_map_add_wrap;
3682 sol_map->sol.add_empty = track_empty ? &sol_map_add_empty_wrap : NULL;
3683 space = isl_space_copy(sol_map->sol.space);
3684 sol_map->map = isl_map_alloc_space(space, 1, ISL_MAP_DISJOINT);
3685 if (!sol_map->map)
3686 goto error;
3688 if (track_empty) {
3689 sol_map->empty = isl_set_alloc_space(isl_basic_set_get_space(dom),
3690 1, ISL_SET_DISJOINT);
3691 if (!sol_map->empty)
3692 goto error;
3695 isl_basic_set_free(dom);
3696 return &sol_map->sol;
3697 error:
3698 isl_basic_set_free(dom);
3699 sol_free(&sol_map->sol);
3700 return NULL;
3703 /* Check whether all coefficients of (non-parameter) variables
3704 * are non-positive, meaning that no pivots can be performed on the row.
3706 static int is_critical(struct isl_tab *tab, int row)
3708 int j;
3709 unsigned off = 2 + tab->M;
3711 for (j = tab->n_dead; j < tab->n_col; ++j) {
3712 if (col_is_parameter_var(tab, j))
3713 continue;
3715 if (isl_int_is_pos(tab->mat->row[row][off + j]))
3716 return 0;
3719 return 1;
3722 /* Check whether the inequality represented by vec is strict over the integers,
3723 * i.e., there are no integer values satisfying the constraint with
3724 * equality. This happens if the gcd of the coefficients is not a divisor
3725 * of the constant term. If so, scale the constraint down by the gcd
3726 * of the coefficients.
3728 static int is_strict(struct isl_vec *vec)
3730 isl_int gcd;
3731 int strict = 0;
3733 isl_int_init(gcd);
3734 isl_seq_gcd(vec->el + 1, vec->size - 1, &gcd);
3735 if (!isl_int_is_one(gcd)) {
3736 strict = !isl_int_is_divisible_by(vec->el[0], gcd);
3737 isl_int_fdiv_q(vec->el[0], vec->el[0], gcd);
3738 isl_seq_scale_down(vec->el + 1, vec->el + 1, gcd, vec->size-1);
3740 isl_int_clear(gcd);
3742 return strict;
3745 /* Determine the sign of the given row of the main tableau.
3746 * The result is one of
3747 * isl_tab_row_pos: always non-negative; no pivot needed
3748 * isl_tab_row_neg: always non-positive; pivot
3749 * isl_tab_row_any: can be both positive and negative; split
3751 * We first handle some simple cases
3752 * - the row sign may be known already
3753 * - the row may be obviously non-negative
3754 * - the parametric constant may be equal to that of another row
3755 * for which we know the sign. This sign will be either "pos" or
3756 * "any". If it had been "neg" then we would have pivoted before.
3758 * If none of these cases hold, we check the value of the row for each
3759 * of the currently active samples. Based on the signs of these values
3760 * we make an initial determination of the sign of the row.
3762 * all zero -> unk(nown)
3763 * all non-negative -> pos
3764 * all non-positive -> neg
3765 * both negative and positive -> all
3767 * If we end up with "all", we are done.
3768 * Otherwise, we perform a check for positive and/or negative
3769 * values as follows.
3771 * samples neg unk pos
3772 * <0 ? Y N Y N
3773 * pos any pos
3774 * >0 ? Y N Y N
3775 * any neg any neg
3777 * There is no special sign for "zero", because we can usually treat zero
3778 * as either non-negative or non-positive, whatever works out best.
3779 * However, if the row is "critical", meaning that pivoting is impossible
3780 * then we don't want to limp zero with the non-positive case, because
3781 * then we we would lose the solution for those values of the parameters
3782 * where the value of the row is zero. Instead, we treat 0 as non-negative
3783 * ensuring a split if the row can attain both zero and negative values.
3784 * The same happens when the original constraint was one that could not
3785 * be satisfied with equality by any integer values of the parameters.
3786 * In this case, we normalize the constraint, but then a value of zero
3787 * for the normalized constraint is actually a positive value for the
3788 * original constraint, so again we need to treat zero as non-negative.
3789 * In both these cases, we have the following decision tree instead:
3791 * all non-negative -> pos
3792 * all negative -> neg
3793 * both negative and non-negative -> all
3795 * samples neg pos
3796 * <0 ? Y N
3797 * any pos
3798 * >=0 ? Y N
3799 * any neg
3801 static enum isl_tab_row_sign row_sign(struct isl_tab *tab,
3802 struct isl_sol *sol, int row)
3804 struct isl_vec *ineq = NULL;
3805 enum isl_tab_row_sign res = isl_tab_row_unknown;
3806 int critical;
3807 int strict;
3808 int row2;
3810 if (tab->row_sign[row] != isl_tab_row_unknown)
3811 return tab->row_sign[row];
3812 if (is_obviously_nonneg(tab, row))
3813 return isl_tab_row_pos;
3814 for (row2 = tab->n_redundant; row2 < tab->n_row; ++row2) {
3815 if (tab->row_sign[row2] == isl_tab_row_unknown)
3816 continue;
3817 if (identical_parameter_line(tab, row, row2))
3818 return tab->row_sign[row2];
3821 critical = is_critical(tab, row);
3823 ineq = get_row_parameter_ineq(tab, row);
3824 if (!ineq)
3825 goto error;
3827 strict = is_strict(ineq);
3829 res = sol->context->op->ineq_sign(sol->context, ineq->el,
3830 critical || strict);
3832 if (res == isl_tab_row_unknown || res == isl_tab_row_pos) {
3833 /* test for negative values */
3834 int feasible;
3835 isl_seq_neg(ineq->el, ineq->el, ineq->size);
3836 isl_int_sub_ui(ineq->el[0], ineq->el[0], 1);
3838 feasible = sol->context->op->test_ineq(sol->context, ineq->el);
3839 if (feasible < 0)
3840 goto error;
3841 if (!feasible)
3842 res = isl_tab_row_pos;
3843 else
3844 res = (res == isl_tab_row_unknown) ? isl_tab_row_neg
3845 : isl_tab_row_any;
3846 if (res == isl_tab_row_neg) {
3847 isl_seq_neg(ineq->el, ineq->el, ineq->size);
3848 isl_int_sub_ui(ineq->el[0], ineq->el[0], 1);
3852 if (res == isl_tab_row_neg) {
3853 /* test for positive values */
3854 int feasible;
3855 if (!critical && !strict)
3856 isl_int_sub_ui(ineq->el[0], ineq->el[0], 1);
3858 feasible = sol->context->op->test_ineq(sol->context, ineq->el);
3859 if (feasible < 0)
3860 goto error;
3861 if (feasible)
3862 res = isl_tab_row_any;
3865 isl_vec_free(ineq);
3866 return res;
3867 error:
3868 isl_vec_free(ineq);
3869 return isl_tab_row_unknown;
3872 static void find_solutions(struct isl_sol *sol, struct isl_tab *tab);
3874 /* Find solutions for values of the parameters that satisfy the given
3875 * inequality.
3877 * We currently take a snapshot of the context tableau that is reset
3878 * when we return from this function, while we make a copy of the main
3879 * tableau, leaving the original main tableau untouched.
3880 * These are fairly arbitrary choices. Making a copy also of the context
3881 * tableau would obviate the need to undo any changes made to it later,
3882 * while taking a snapshot of the main tableau could reduce memory usage.
3883 * If we were to switch to taking a snapshot of the main tableau,
3884 * we would have to keep in mind that we need to save the row signs
3885 * and that we need to do this before saving the current basis
3886 * such that the basis has been restore before we restore the row signs.
3888 static void find_in_pos(struct isl_sol *sol, struct isl_tab *tab, isl_int *ineq)
3890 void *saved;
3892 if (!sol->context)
3893 goto error;
3894 saved = sol->context->op->save(sol->context);
3896 tab = isl_tab_dup(tab);
3897 if (!tab)
3898 goto error;
3900 sol->context->op->add_ineq(sol->context, ineq, 0, 1);
3902 find_solutions(sol, tab);
3904 if (!sol->error)
3905 sol->context->op->restore(sol->context, saved);
3906 else
3907 sol->context->op->discard(saved);
3908 return;
3909 error:
3910 sol->error = 1;
3913 /* Record the absence of solutions for those values of the parameters
3914 * that do not satisfy the given inequality with equality.
3916 static void no_sol_in_strict(struct isl_sol *sol,
3917 struct isl_tab *tab, struct isl_vec *ineq)
3919 int empty;
3920 void *saved;
3922 if (!sol->context || sol->error)
3923 goto error;
3924 saved = sol->context->op->save(sol->context);
3926 isl_int_sub_ui(ineq->el[0], ineq->el[0], 1);
3928 sol->context->op->add_ineq(sol->context, ineq->el, 1, 0);
3929 if (!sol->context)
3930 goto error;
3932 empty = tab->empty;
3933 tab->empty = 1;
3934 sol_add(sol, tab);
3935 tab->empty = empty;
3937 isl_int_add_ui(ineq->el[0], ineq->el[0], 1);
3939 sol->context->op->restore(sol->context, saved);
3940 return;
3941 error:
3942 sol->error = 1;
3945 /* Reset all row variables that are marked to have a sign that may
3946 * be both positive and negative to have an unknown sign.
3948 static void reset_any_to_unknown(struct isl_tab *tab)
3950 int row;
3952 for (row = tab->n_redundant; row < tab->n_row; ++row) {
3953 if (!isl_tab_var_from_row(tab, row)->is_nonneg)
3954 continue;
3955 if (tab->row_sign[row] == isl_tab_row_any)
3956 tab->row_sign[row] = isl_tab_row_unknown;
3960 /* Compute the lexicographic minimum of the set represented by the main
3961 * tableau "tab" within the context "sol->context_tab".
3962 * On entry the sample value of the main tableau is lexicographically
3963 * less than or equal to this lexicographic minimum.
3964 * Pivots are performed until a feasible point is found, which is then
3965 * necessarily equal to the minimum, or until the tableau is found to
3966 * be infeasible. Some pivots may need to be performed for only some
3967 * feasible values of the context tableau. If so, the context tableau
3968 * is split into a part where the pivot is needed and a part where it is not.
3970 * Whenever we enter the main loop, the main tableau is such that no
3971 * "obvious" pivots need to be performed on it, where "obvious" means
3972 * that the given row can be seen to be negative without looking at
3973 * the context tableau. In particular, for non-parametric problems,
3974 * no pivots need to be performed on the main tableau.
3975 * The caller of find_solutions is responsible for making this property
3976 * hold prior to the first iteration of the loop, while restore_lexmin
3977 * is called before every other iteration.
3979 * Inside the main loop, we first examine the signs of the rows of
3980 * the main tableau within the context of the context tableau.
3981 * If we find a row that is always non-positive for all values of
3982 * the parameters satisfying the context tableau and negative for at
3983 * least one value of the parameters, we perform the appropriate pivot
3984 * and start over. An exception is the case where no pivot can be
3985 * performed on the row. In this case, we require that the sign of
3986 * the row is negative for all values of the parameters (rather than just
3987 * non-positive). This special case is handled inside row_sign, which
3988 * will say that the row can have any sign if it determines that it can
3989 * attain both negative and zero values.
3991 * If we can't find a row that always requires a pivot, but we can find
3992 * one or more rows that require a pivot for some values of the parameters
3993 * (i.e., the row can attain both positive and negative signs), then we split
3994 * the context tableau into two parts, one where we force the sign to be
3995 * non-negative and one where we force is to be negative.
3996 * The non-negative part is handled by a recursive call (through find_in_pos).
3997 * Upon returning from this call, we continue with the negative part and
3998 * perform the required pivot.
4000 * If no such rows can be found, all rows are non-negative and we have
4001 * found a (rational) feasible point. If we only wanted a rational point
4002 * then we are done.
4003 * Otherwise, we check if all values of the sample point of the tableau
4004 * are integral for the variables. If so, we have found the minimal
4005 * integral point and we are done.
4006 * If the sample point is not integral, then we need to make a distinction
4007 * based on whether the constant term is non-integral or the coefficients
4008 * of the parameters. Furthermore, in order to decide how to handle
4009 * the non-integrality, we also need to know whether the coefficients
4010 * of the other columns in the tableau are integral. This leads
4011 * to the following table. The first two rows do not correspond
4012 * to a non-integral sample point and are only mentioned for completeness.
4014 * constant parameters other
4016 * int int int |
4017 * int int rat | -> no problem
4019 * rat int int -> fail
4021 * rat int rat -> cut
4023 * int rat rat |
4024 * rat rat rat | -> parametric cut
4026 * int rat int |
4027 * rat rat int | -> split context
4029 * If the parametric constant is completely integral, then there is nothing
4030 * to be done. If the constant term is non-integral, but all the other
4031 * coefficient are integral, then there is nothing that can be done
4032 * and the tableau has no integral solution.
4033 * If, on the other hand, one or more of the other columns have rational
4034 * coefficients, but the parameter coefficients are all integral, then
4035 * we can perform a regular (non-parametric) cut.
4036 * Finally, if there is any parameter coefficient that is non-integral,
4037 * then we need to involve the context tableau. There are two cases here.
4038 * If at least one other column has a rational coefficient, then we
4039 * can perform a parametric cut in the main tableau by adding a new
4040 * integer division in the context tableau.
4041 * If all other columns have integral coefficients, then we need to
4042 * enforce that the rational combination of parameters (c + \sum a_i y_i)/m
4043 * is always integral. We do this by introducing an integer division
4044 * q = floor((c + \sum a_i y_i)/m) and stipulating that its argument should
4045 * always be integral in the context tableau, i.e., m q = c + \sum a_i y_i.
4046 * Since q is expressed in the tableau as
4047 * c + \sum a_i y_i - m q >= 0
4048 * -c - \sum a_i y_i + m q + m - 1 >= 0
4049 * it is sufficient to add the inequality
4050 * -c - \sum a_i y_i + m q >= 0
4051 * In the part of the context where this inequality does not hold, the
4052 * main tableau is marked as being empty.
4054 static void find_solutions(struct isl_sol *sol, struct isl_tab *tab)
4056 struct isl_context *context;
4057 int r;
4059 if (!tab || sol->error)
4060 goto error;
4062 context = sol->context;
4064 if (tab->empty)
4065 goto done;
4066 if (context->op->is_empty(context))
4067 goto done;
4069 for (r = 0; r >= 0 && tab && !tab->empty; r = restore_lexmin(tab)) {
4070 int flags;
4071 int row;
4072 enum isl_tab_row_sign sgn;
4073 int split = -1;
4074 int n_split = 0;
4076 for (row = tab->n_redundant; row < tab->n_row; ++row) {
4077 if (!isl_tab_var_from_row(tab, row)->is_nonneg)
4078 continue;
4079 sgn = row_sign(tab, sol, row);
4080 if (!sgn)
4081 goto error;
4082 tab->row_sign[row] = sgn;
4083 if (sgn == isl_tab_row_any)
4084 n_split++;
4085 if (sgn == isl_tab_row_any && split == -1)
4086 split = row;
4087 if (sgn == isl_tab_row_neg)
4088 break;
4090 if (row < tab->n_row)
4091 continue;
4092 if (split != -1) {
4093 struct isl_vec *ineq;
4094 if (n_split != 1)
4095 split = context->op->best_split(context, tab);
4096 if (split < 0)
4097 goto error;
4098 ineq = get_row_parameter_ineq(tab, split);
4099 if (!ineq)
4100 goto error;
4101 is_strict(ineq);
4102 reset_any_to_unknown(tab);
4103 tab->row_sign[split] = isl_tab_row_pos;
4104 sol_inc_level(sol);
4105 find_in_pos(sol, tab, ineq->el);
4106 tab->row_sign[split] = isl_tab_row_neg;
4107 isl_seq_neg(ineq->el, ineq->el, ineq->size);
4108 isl_int_sub_ui(ineq->el[0], ineq->el[0], 1);
4109 if (!sol->error)
4110 context->op->add_ineq(context, ineq->el, 0, 1);
4111 isl_vec_free(ineq);
4112 if (sol->error)
4113 goto error;
4114 continue;
4116 if (tab->rational)
4117 break;
4118 row = first_non_integer_row(tab, &flags);
4119 if (row < 0)
4120 break;
4121 if (ISL_FL_ISSET(flags, I_PAR)) {
4122 if (ISL_FL_ISSET(flags, I_VAR)) {
4123 if (isl_tab_mark_empty(tab) < 0)
4124 goto error;
4125 break;
4127 row = add_cut(tab, row);
4128 } else if (ISL_FL_ISSET(flags, I_VAR)) {
4129 struct isl_vec *div;
4130 struct isl_vec *ineq;
4131 int d;
4132 div = get_row_split_div(tab, row);
4133 if (!div)
4134 goto error;
4135 d = context->op->get_div(context, tab, div);
4136 isl_vec_free(div);
4137 if (d < 0)
4138 goto error;
4139 ineq = ineq_for_div(context->op->peek_basic_set(context), d);
4140 if (!ineq)
4141 goto error;
4142 sol_inc_level(sol);
4143 no_sol_in_strict(sol, tab, ineq);
4144 isl_seq_neg(ineq->el, ineq->el, ineq->size);
4145 context->op->add_ineq(context, ineq->el, 1, 1);
4146 isl_vec_free(ineq);
4147 if (sol->error || !context->op->is_ok(context))
4148 goto error;
4149 tab = set_row_cst_to_div(tab, row, d);
4150 if (context->op->is_empty(context))
4151 break;
4152 } else
4153 row = add_parametric_cut(tab, row, context);
4154 if (row < 0)
4155 goto error;
4157 if (r < 0)
4158 goto error;
4159 done:
4160 sol_add(sol, tab);
4161 isl_tab_free(tab);
4162 return;
4163 error:
4164 isl_tab_free(tab);
4165 sol->error = 1;
4168 /* Does "sol" contain a pair of partial solutions that could potentially
4169 * be merged?
4171 * We currently only check that "sol" is not in an error state
4172 * and that there are at least two partial solutions of which the final two
4173 * are defined at the same level.
4175 static int sol_has_mergeable_solutions(struct isl_sol *sol)
4177 if (sol->error)
4178 return 0;
4179 if (!sol->partial)
4180 return 0;
4181 if (!sol->partial->next)
4182 return 0;
4183 return sol->partial->level == sol->partial->next->level;
4186 /* Compute the lexicographic minimum of the set represented by the main
4187 * tableau "tab" within the context "sol->context_tab".
4189 * As a preprocessing step, we first transfer all the purely parametric
4190 * equalities from the main tableau to the context tableau, i.e.,
4191 * parameters that have been pivoted to a row.
4192 * These equalities are ignored by the main algorithm, because the
4193 * corresponding rows may not be marked as being non-negative.
4194 * In parts of the context where the added equality does not hold,
4195 * the main tableau is marked as being empty.
4197 * Before we embark on the actual computation, we save a copy
4198 * of the context. When we return, we check if there are any
4199 * partial solutions that can potentially be merged. If so,
4200 * we perform a rollback to the initial state of the context.
4201 * The merging of partial solutions happens inside calls to
4202 * sol_dec_level that are pushed onto the undo stack of the context.
4203 * If there are no partial solutions that can potentially be merged
4204 * then the rollback is skipped as it would just be wasted effort.
4206 static void find_solutions_main(struct isl_sol *sol, struct isl_tab *tab)
4208 int row;
4209 void *saved;
4211 if (!tab)
4212 goto error;
4214 sol->level = 0;
4216 for (row = tab->n_redundant; row < tab->n_row; ++row) {
4217 int p;
4218 struct isl_vec *eq;
4220 if (!row_is_parameter_var(tab, row))
4221 continue;
4222 if (tab->row_var[row] < tab->n_param)
4223 p = tab->row_var[row];
4224 else
4225 p = tab->row_var[row]
4226 + tab->n_param - (tab->n_var - tab->n_div);
4228 eq = isl_vec_alloc(tab->mat->ctx, 1+tab->n_param+tab->n_div);
4229 if (!eq)
4230 goto error;
4231 get_row_parameter_line(tab, row, eq->el);
4232 isl_int_neg(eq->el[1 + p], tab->mat->row[row][0]);
4233 eq = isl_vec_normalize(eq);
4235 sol_inc_level(sol);
4236 no_sol_in_strict(sol, tab, eq);
4238 isl_seq_neg(eq->el, eq->el, eq->size);
4239 sol_inc_level(sol);
4240 no_sol_in_strict(sol, tab, eq);
4241 isl_seq_neg(eq->el, eq->el, eq->size);
4243 sol->context->op->add_eq(sol->context, eq->el, 1, 1);
4245 isl_vec_free(eq);
4247 if (isl_tab_mark_redundant(tab, row) < 0)
4248 goto error;
4250 if (sol->context->op->is_empty(sol->context))
4251 break;
4253 row = tab->n_redundant - 1;
4256 saved = sol->context->op->save(sol->context);
4258 find_solutions(sol, tab);
4260 if (sol_has_mergeable_solutions(sol))
4261 sol->context->op->restore(sol->context, saved);
4262 else
4263 sol->context->op->discard(saved);
4265 sol->level = 0;
4266 sol_pop(sol);
4268 return;
4269 error:
4270 isl_tab_free(tab);
4271 sol->error = 1;
4274 /* Check if integer division "div" of "dom" also occurs in "bmap".
4275 * If so, return its position within the divs.
4276 * If not, return -1.
4278 static int find_context_div(struct isl_basic_map *bmap,
4279 struct isl_basic_set *dom, unsigned div)
4281 int i;
4282 unsigned b_dim = isl_space_dim(bmap->dim, isl_dim_all);
4283 unsigned d_dim = isl_space_dim(dom->dim, isl_dim_all);
4285 if (isl_int_is_zero(dom->div[div][0]))
4286 return -1;
4287 if (isl_seq_first_non_zero(dom->div[div] + 2 + d_dim, dom->n_div) != -1)
4288 return -1;
4290 for (i = 0; i < bmap->n_div; ++i) {
4291 if (isl_int_is_zero(bmap->div[i][0]))
4292 continue;
4293 if (isl_seq_first_non_zero(bmap->div[i] + 2 + d_dim,
4294 (b_dim - d_dim) + bmap->n_div) != -1)
4295 continue;
4296 if (isl_seq_eq(bmap->div[i], dom->div[div], 2 + d_dim))
4297 return i;
4299 return -1;
4302 /* The correspondence between the variables in the main tableau,
4303 * the context tableau, and the input map and domain is as follows.
4304 * The first n_param and the last n_div variables of the main tableau
4305 * form the variables of the context tableau.
4306 * In the basic map, these n_param variables correspond to the
4307 * parameters and the input dimensions. In the domain, they correspond
4308 * to the parameters and the set dimensions.
4309 * The n_div variables correspond to the integer divisions in the domain.
4310 * To ensure that everything lines up, we may need to copy some of the
4311 * integer divisions of the domain to the map. These have to be placed
4312 * in the same order as those in the context and they have to be placed
4313 * after any other integer divisions that the map may have.
4314 * This function performs the required reordering.
4316 static __isl_give isl_basic_map *align_context_divs(
4317 __isl_take isl_basic_map *bmap, __isl_keep isl_basic_set *dom)
4319 int i;
4320 int common = 0;
4321 int other;
4323 for (i = 0; i < dom->n_div; ++i)
4324 if (find_context_div(bmap, dom, i) != -1)
4325 common++;
4326 other = bmap->n_div - common;
4327 if (dom->n_div - common > 0) {
4328 bmap = isl_basic_map_extend_space(bmap, isl_space_copy(bmap->dim),
4329 dom->n_div - common, 0, 0);
4330 if (!bmap)
4331 return NULL;
4333 for (i = 0; i < dom->n_div; ++i) {
4334 int pos = find_context_div(bmap, dom, i);
4335 if (pos < 0) {
4336 pos = isl_basic_map_alloc_div(bmap);
4337 if (pos < 0)
4338 goto error;
4339 isl_int_set_si(bmap->div[pos][0], 0);
4341 if (pos != other + i)
4342 isl_basic_map_swap_div(bmap, pos, other + i);
4344 return bmap;
4345 error:
4346 isl_basic_map_free(bmap);
4347 return NULL;
4350 /* Base case of isl_tab_basic_map_partial_lexopt, after removing
4351 * some obvious symmetries.
4353 * We make sure the divs in the domain are properly ordered,
4354 * because they will be added one by one in the given order
4355 * during the construction of the solution map.
4356 * Furthermore, make sure that the known integer divisions
4357 * appear before any unknown integer division because the solution
4358 * may depend on the known integer divisions, while anything that
4359 * depends on any variable starting from the first unknown integer
4360 * division is ignored in sol_pma_add.
4362 static struct isl_sol *basic_map_partial_lexopt_base_sol(
4363 __isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom,
4364 __isl_give isl_set **empty, int max,
4365 struct isl_sol *(*init)(__isl_keep isl_basic_map *bmap,
4366 __isl_take isl_basic_set *dom, int track_empty, int max))
4368 struct isl_tab *tab;
4369 struct isl_sol *sol = NULL;
4370 struct isl_context *context;
4372 if (dom->n_div) {
4373 dom = isl_basic_set_sort_divs(dom);
4374 bmap = align_context_divs(bmap, dom);
4376 sol = init(bmap, dom, !!empty, max);
4377 if (!sol)
4378 goto error;
4380 context = sol->context;
4381 if (isl_basic_set_plain_is_empty(context->op->peek_basic_set(context)))
4382 /* nothing */;
4383 else if (isl_basic_map_plain_is_empty(bmap)) {
4384 if (sol->add_empty)
4385 sol->add_empty(sol,
4386 isl_basic_set_copy(context->op->peek_basic_set(context)));
4387 } else {
4388 tab = tab_for_lexmin(bmap,
4389 context->op->peek_basic_set(context), 1, max);
4390 tab = context->op->detect_nonnegative_parameters(context, tab);
4391 find_solutions_main(sol, tab);
4393 if (sol->error)
4394 goto error;
4396 isl_basic_map_free(bmap);
4397 return sol;
4398 error:
4399 sol_free(sol);
4400 isl_basic_map_free(bmap);
4401 return NULL;
4404 /* Base case of isl_tab_basic_map_partial_lexopt, after removing
4405 * some obvious symmetries.
4407 * We call basic_map_partial_lexopt_base_sol and extract the results.
4409 static __isl_give isl_map *basic_map_partial_lexopt_base(
4410 __isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom,
4411 __isl_give isl_set **empty, int max)
4413 isl_map *result = NULL;
4414 struct isl_sol *sol;
4415 struct isl_sol_map *sol_map;
4417 sol = basic_map_partial_lexopt_base_sol(bmap, dom, empty, max,
4418 &sol_map_init);
4419 if (!sol)
4420 return NULL;
4421 sol_map = (struct isl_sol_map *) sol;
4423 result = isl_map_copy(sol_map->map);
4424 if (empty)
4425 *empty = isl_set_copy(sol_map->empty);
4426 sol_free(&sol_map->sol);
4427 return result;
4430 /* Return a count of the number of occurrences of the "n" first
4431 * variables in the inequality constraints of "bmap".
4433 static __isl_give int *count_occurrences(__isl_keep isl_basic_map *bmap,
4434 int n)
4436 int i, j;
4437 isl_ctx *ctx;
4438 int *occurrences;
4440 if (!bmap)
4441 return NULL;
4442 ctx = isl_basic_map_get_ctx(bmap);
4443 occurrences = isl_calloc_array(ctx, int, n);
4444 if (!occurrences)
4445 return NULL;
4447 for (i = 0; i < bmap->n_ineq; ++i) {
4448 for (j = 0; j < n; ++j) {
4449 if (!isl_int_is_zero(bmap->ineq[i][1 + j]))
4450 occurrences[j]++;
4454 return occurrences;
4457 /* Do all of the "n" variables with non-zero coefficients in "c"
4458 * occur in exactly a single constraint.
4459 * "occurrences" is an array of length "n" containing the number
4460 * of occurrences of each of the variables in the inequality constraints.
4462 static int single_occurrence(int n, isl_int *c, int *occurrences)
4464 int i;
4466 for (i = 0; i < n; ++i) {
4467 if (isl_int_is_zero(c[i]))
4468 continue;
4469 if (occurrences[i] != 1)
4470 return 0;
4473 return 1;
4476 /* Do all of the "n" initial variables that occur in inequality constraint
4477 * "ineq" of "bmap" only occur in that constraint?
4479 static int all_single_occurrence(__isl_keep isl_basic_map *bmap, int ineq,
4480 int n)
4482 int i, j;
4484 for (i = 0; i < n; ++i) {
4485 if (isl_int_is_zero(bmap->ineq[ineq][1 + i]))
4486 continue;
4487 for (j = 0; j < bmap->n_ineq; ++j) {
4488 if (j == ineq)
4489 continue;
4490 if (!isl_int_is_zero(bmap->ineq[j][1 + i]))
4491 return 0;
4495 return 1;
4498 /* Structure used during detection of parallel constraints.
4499 * n_in: number of "input" variables: isl_dim_param + isl_dim_in
4500 * n_out: number of "output" variables: isl_dim_out + isl_dim_div
4501 * val: the coefficients of the output variables
4503 struct isl_constraint_equal_info {
4504 unsigned n_in;
4505 unsigned n_out;
4506 isl_int *val;
4509 /* Check whether the coefficients of the output variables
4510 * of the constraint in "entry" are equal to info->val.
4512 static int constraint_equal(const void *entry, const void *val)
4514 isl_int **row = (isl_int **)entry;
4515 const struct isl_constraint_equal_info *info = val;
4517 return isl_seq_eq((*row) + 1 + info->n_in, info->val, info->n_out);
4520 /* Check whether "bmap" has a pair of constraints that have
4521 * the same coefficients for the output variables.
4522 * Note that the coefficients of the existentially quantified
4523 * variables need to be zero since the existentially quantified
4524 * of the result are usually not the same as those of the input.
4525 * Furthermore, check that each of the input variables that occur
4526 * in those constraints does not occur in any other constraint.
4527 * If so, return true and return the row indices of the two constraints
4528 * in *first and *second.
4530 static isl_bool parallel_constraints(__isl_keep isl_basic_map *bmap,
4531 int *first, int *second)
4533 int i;
4534 isl_ctx *ctx;
4535 int *occurrences = NULL;
4536 struct isl_hash_table *table = NULL;
4537 struct isl_hash_table_entry *entry;
4538 struct isl_constraint_equal_info info;
4539 unsigned n_out;
4540 unsigned n_div;
4542 ctx = isl_basic_map_get_ctx(bmap);
4543 table = isl_hash_table_alloc(ctx, bmap->n_ineq);
4544 if (!table)
4545 goto error;
4547 info.n_in = isl_basic_map_dim(bmap, isl_dim_param) +
4548 isl_basic_map_dim(bmap, isl_dim_in);
4549 occurrences = count_occurrences(bmap, info.n_in);
4550 if (info.n_in && !occurrences)
4551 goto error;
4552 n_out = isl_basic_map_dim(bmap, isl_dim_out);
4553 n_div = isl_basic_map_dim(bmap, isl_dim_div);
4554 info.n_out = n_out + n_div;
4555 for (i = 0; i < bmap->n_ineq; ++i) {
4556 uint32_t hash;
4558 info.val = bmap->ineq[i] + 1 + info.n_in;
4559 if (isl_seq_first_non_zero(info.val, n_out) < 0)
4560 continue;
4561 if (isl_seq_first_non_zero(info.val + n_out, n_div) >= 0)
4562 continue;
4563 if (!single_occurrence(info.n_in, bmap->ineq[i] + 1,
4564 occurrences))
4565 continue;
4566 hash = isl_seq_get_hash(info.val, info.n_out);
4567 entry = isl_hash_table_find(ctx, table, hash,
4568 constraint_equal, &info, 1);
4569 if (!entry)
4570 goto error;
4571 if (entry->data)
4572 break;
4573 entry->data = &bmap->ineq[i];
4576 if (i < bmap->n_ineq) {
4577 *first = ((isl_int **)entry->data) - bmap->ineq;
4578 *second = i;
4581 isl_hash_table_free(ctx, table);
4582 free(occurrences);
4584 return i < bmap->n_ineq;
4585 error:
4586 isl_hash_table_free(ctx, table);
4587 free(occurrences);
4588 return isl_bool_error;
4591 /* Given a set of upper bounds in "var", add constraints to "bset"
4592 * that make the i-th bound smallest.
4594 * In particular, if there are n bounds b_i, then add the constraints
4596 * b_i <= b_j for j > i
4597 * b_i < b_j for j < i
4599 static __isl_give isl_basic_set *select_minimum(__isl_take isl_basic_set *bset,
4600 __isl_keep isl_mat *var, int i)
4602 isl_ctx *ctx;
4603 int j, k;
4605 ctx = isl_mat_get_ctx(var);
4607 for (j = 0; j < var->n_row; ++j) {
4608 if (j == i)
4609 continue;
4610 k = isl_basic_set_alloc_inequality(bset);
4611 if (k < 0)
4612 goto error;
4613 isl_seq_combine(bset->ineq[k], ctx->one, var->row[j],
4614 ctx->negone, var->row[i], var->n_col);
4615 isl_int_set_si(bset->ineq[k][var->n_col], 0);
4616 if (j < i)
4617 isl_int_sub_ui(bset->ineq[k][0], bset->ineq[k][0], 1);
4620 bset = isl_basic_set_finalize(bset);
4622 return bset;
4623 error:
4624 isl_basic_set_free(bset);
4625 return NULL;
4628 /* Given a set of upper bounds on the last "input" variable m,
4629 * construct a set that assigns the minimal upper bound to m, i.e.,
4630 * construct a set that divides the space into cells where one
4631 * of the upper bounds is smaller than all the others and assign
4632 * this upper bound to m.
4634 * In particular, if there are n bounds b_i, then the result
4635 * consists of n basic sets, each one of the form
4637 * m = b_i
4638 * b_i <= b_j for j > i
4639 * b_i < b_j for j < i
4641 static __isl_give isl_set *set_minimum(__isl_take isl_space *dim,
4642 __isl_take isl_mat *var)
4644 int i, k;
4645 isl_basic_set *bset = NULL;
4646 isl_set *set = NULL;
4648 if (!dim || !var)
4649 goto error;
4651 set = isl_set_alloc_space(isl_space_copy(dim),
4652 var->n_row, ISL_SET_DISJOINT);
4654 for (i = 0; i < var->n_row; ++i) {
4655 bset = isl_basic_set_alloc_space(isl_space_copy(dim), 0,
4656 1, var->n_row - 1);
4657 k = isl_basic_set_alloc_equality(bset);
4658 if (k < 0)
4659 goto error;
4660 isl_seq_cpy(bset->eq[k], var->row[i], var->n_col);
4661 isl_int_set_si(bset->eq[k][var->n_col], -1);
4662 bset = select_minimum(bset, var, i);
4663 set = isl_set_add_basic_set(set, bset);
4666 isl_space_free(dim);
4667 isl_mat_free(var);
4668 return set;
4669 error:
4670 isl_basic_set_free(bset);
4671 isl_set_free(set);
4672 isl_space_free(dim);
4673 isl_mat_free(var);
4674 return NULL;
4677 /* Given that the last input variable of "bmap" represents the minimum
4678 * of the bounds in "cst", check whether we need to split the domain
4679 * based on which bound attains the minimum.
4681 * A split is needed when the minimum appears in an integer division
4682 * or in an equality. Otherwise, it is only needed if it appears in
4683 * an upper bound that is different from the upper bounds on which it
4684 * is defined.
4686 static isl_bool need_split_basic_map(__isl_keep isl_basic_map *bmap,
4687 __isl_keep isl_mat *cst)
4689 int i, j;
4690 unsigned total;
4691 unsigned pos;
4693 pos = cst->n_col - 1;
4694 total = isl_basic_map_dim(bmap, isl_dim_all);
4696 for (i = 0; i < bmap->n_div; ++i)
4697 if (!isl_int_is_zero(bmap->div[i][2 + pos]))
4698 return isl_bool_true;
4700 for (i = 0; i < bmap->n_eq; ++i)
4701 if (!isl_int_is_zero(bmap->eq[i][1 + pos]))
4702 return isl_bool_true;
4704 for (i = 0; i < bmap->n_ineq; ++i) {
4705 if (isl_int_is_nonneg(bmap->ineq[i][1 + pos]))
4706 continue;
4707 if (!isl_int_is_negone(bmap->ineq[i][1 + pos]))
4708 return isl_bool_true;
4709 if (isl_seq_first_non_zero(bmap->ineq[i] + 1 + pos + 1,
4710 total - pos - 1) >= 0)
4711 return isl_bool_true;
4713 for (j = 0; j < cst->n_row; ++j)
4714 if (isl_seq_eq(bmap->ineq[i], cst->row[j], cst->n_col))
4715 break;
4716 if (j >= cst->n_row)
4717 return isl_bool_true;
4720 return isl_bool_false;
4723 /* Given that the last set variable of "bset" represents the minimum
4724 * of the bounds in "cst", check whether we need to split the domain
4725 * based on which bound attains the minimum.
4727 * We simply call need_split_basic_map here. This is safe because
4728 * the position of the minimum is computed from "cst" and not
4729 * from "bmap".
4731 static isl_bool need_split_basic_set(__isl_keep isl_basic_set *bset,
4732 __isl_keep isl_mat *cst)
4734 return need_split_basic_map(bset_to_bmap(bset), cst);
4737 /* Given that the last set variable of "set" represents the minimum
4738 * of the bounds in "cst", check whether we need to split the domain
4739 * based on which bound attains the minimum.
4741 static isl_bool need_split_set(__isl_keep isl_set *set, __isl_keep isl_mat *cst)
4743 int i;
4745 for (i = 0; i < set->n; ++i) {
4746 isl_bool split;
4748 split = need_split_basic_set(set->p[i], cst);
4749 if (split < 0 || split)
4750 return split;
4753 return isl_bool_false;
4756 /* Given a set of which the last set variable is the minimum
4757 * of the bounds in "cst", split each basic set in the set
4758 * in pieces where one of the bounds is (strictly) smaller than the others.
4759 * This subdivision is given in "min_expr".
4760 * The variable is subsequently projected out.
4762 * We only do the split when it is needed.
4763 * For example if the last input variable m = min(a,b) and the only
4764 * constraints in the given basic set are lower bounds on m,
4765 * i.e., l <= m = min(a,b), then we can simply project out m
4766 * to obtain l <= a and l <= b, without having to split on whether
4767 * m is equal to a or b.
4769 static __isl_give isl_set *split(__isl_take isl_set *empty,
4770 __isl_take isl_set *min_expr, __isl_take isl_mat *cst)
4772 int n_in;
4773 int i;
4774 isl_space *dim;
4775 isl_set *res;
4777 if (!empty || !min_expr || !cst)
4778 goto error;
4780 n_in = isl_set_dim(empty, isl_dim_set);
4781 dim = isl_set_get_space(empty);
4782 dim = isl_space_drop_dims(dim, isl_dim_set, n_in - 1, 1);
4783 res = isl_set_empty(dim);
4785 for (i = 0; i < empty->n; ++i) {
4786 isl_bool split;
4787 isl_set *set;
4789 set = isl_set_from_basic_set(isl_basic_set_copy(empty->p[i]));
4790 split = need_split_basic_set(empty->p[i], cst);
4791 if (split < 0)
4792 set = isl_set_free(set);
4793 else if (split)
4794 set = isl_set_intersect(set, isl_set_copy(min_expr));
4795 set = isl_set_remove_dims(set, isl_dim_set, n_in - 1, 1);
4797 res = isl_set_union_disjoint(res, set);
4800 isl_set_free(empty);
4801 isl_set_free(min_expr);
4802 isl_mat_free(cst);
4803 return res;
4804 error:
4805 isl_set_free(empty);
4806 isl_set_free(min_expr);
4807 isl_mat_free(cst);
4808 return NULL;
4811 /* Given a map of which the last input variable is the minimum
4812 * of the bounds in "cst", split each basic set in the set
4813 * in pieces where one of the bounds is (strictly) smaller than the others.
4814 * This subdivision is given in "min_expr".
4815 * The variable is subsequently projected out.
4817 * The implementation is essentially the same as that of "split".
4819 static __isl_give isl_map *split_domain(__isl_take isl_map *opt,
4820 __isl_take isl_set *min_expr, __isl_take isl_mat *cst)
4822 int n_in;
4823 int i;
4824 isl_space *dim;
4825 isl_map *res;
4827 if (!opt || !min_expr || !cst)
4828 goto error;
4830 n_in = isl_map_dim(opt, isl_dim_in);
4831 dim = isl_map_get_space(opt);
4832 dim = isl_space_drop_dims(dim, isl_dim_in, n_in - 1, 1);
4833 res = isl_map_empty(dim);
4835 for (i = 0; i < opt->n; ++i) {
4836 isl_map *map;
4837 isl_bool split;
4839 map = isl_map_from_basic_map(isl_basic_map_copy(opt->p[i]));
4840 split = need_split_basic_map(opt->p[i], cst);
4841 if (split < 0)
4842 map = isl_map_free(map);
4843 else if (split)
4844 map = isl_map_intersect_domain(map,
4845 isl_set_copy(min_expr));
4846 map = isl_map_remove_dims(map, isl_dim_in, n_in - 1, 1);
4848 res = isl_map_union_disjoint(res, map);
4851 isl_map_free(opt);
4852 isl_set_free(min_expr);
4853 isl_mat_free(cst);
4854 return res;
4855 error:
4856 isl_map_free(opt);
4857 isl_set_free(min_expr);
4858 isl_mat_free(cst);
4859 return NULL;
4862 static __isl_give isl_map *basic_map_partial_lexopt(
4863 __isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom,
4864 __isl_give isl_set **empty, int max);
4866 /* This function is called from basic_map_partial_lexopt_symm.
4867 * The last variable of "bmap" and "dom" corresponds to the minimum
4868 * of the bounds in "cst". "map_space" is the space of the original
4869 * input relation (of basic_map_partial_lexopt_symm) and "set_space"
4870 * is the space of the original domain.
4872 * We recursively call basic_map_partial_lexopt and then plug in
4873 * the definition of the minimum in the result.
4875 static __isl_give isl_map *basic_map_partial_lexopt_symm_core(
4876 __isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom,
4877 __isl_give isl_set **empty, int max, __isl_take isl_mat *cst,
4878 __isl_take isl_space *map_space, __isl_take isl_space *set_space)
4880 isl_map *opt;
4881 isl_set *min_expr;
4883 min_expr = set_minimum(isl_basic_set_get_space(dom), isl_mat_copy(cst));
4885 opt = basic_map_partial_lexopt(bmap, dom, empty, max);
4887 if (empty) {
4888 *empty = split(*empty,
4889 isl_set_copy(min_expr), isl_mat_copy(cst));
4890 *empty = isl_set_reset_space(*empty, set_space);
4893 opt = split_domain(opt, min_expr, cst);
4894 opt = isl_map_reset_space(opt, map_space);
4896 return opt;
4899 /* Extract a domain from "bmap" for the purpose of computing
4900 * a lexicographic optimum.
4902 * This function is only called when the caller wants to compute a full
4903 * lexicographic optimum, i.e., without specifying a domain. In this case,
4904 * the caller is not interested in the part of the domain space where
4905 * there is no solution and the domain can be initialized to those constraints
4906 * of "bmap" that only involve the parameters and the input dimensions.
4907 * This relieves the parametric programming engine from detecting those
4908 * inequalities and transferring them to the context. More importantly,
4909 * it ensures that those inequalities are transferred first and not
4910 * intermixed with inequalities that actually split the domain.
4912 * If the caller does not require the absence of existentially quantified
4913 * variables in the result (i.e., if ISL_OPT_QE is not set in "flags"),
4914 * then the actual domain of "bmap" can be used. This ensures that
4915 * the domain does not need to be split at all just to separate out
4916 * pieces of the domain that do not have a solution from piece that do.
4917 * This domain cannot be used in general because it may involve
4918 * (unknown) existentially quantified variables which will then also
4919 * appear in the solution.
4921 static __isl_give isl_basic_set *extract_domain(__isl_keep isl_basic_map *bmap,
4922 unsigned flags)
4924 int n_div;
4925 int n_out;
4927 n_div = isl_basic_map_dim(bmap, isl_dim_div);
4928 n_out = isl_basic_map_dim(bmap, isl_dim_out);
4929 bmap = isl_basic_map_copy(bmap);
4930 if (ISL_FL_ISSET(flags, ISL_OPT_QE)) {
4931 bmap = isl_basic_map_drop_constraints_involving_dims(bmap,
4932 isl_dim_div, 0, n_div);
4933 bmap = isl_basic_map_drop_constraints_involving_dims(bmap,
4934 isl_dim_out, 0, n_out);
4936 return isl_basic_map_domain(bmap);
4939 #undef TYPE
4940 #define TYPE isl_map
4941 #undef SUFFIX
4942 #define SUFFIX
4943 #include "isl_tab_lexopt_templ.c"
4945 /* Extract the subsequence of the sample value of "tab"
4946 * starting at "pos" and of length "len".
4948 static __isl_give isl_vec *extract_sample_sequence(struct isl_tab *tab,
4949 int pos, int len)
4951 int i;
4952 isl_ctx *ctx;
4953 isl_vec *v;
4955 ctx = isl_tab_get_ctx(tab);
4956 v = isl_vec_alloc(ctx, len);
4957 if (!v)
4958 return NULL;
4959 for (i = 0; i < len; ++i) {
4960 if (!tab->var[pos + i].is_row) {
4961 isl_int_set_si(v->el[i], 0);
4962 } else {
4963 int row;
4965 row = tab->var[pos + i].index;
4966 isl_int_divexact(v->el[i], tab->mat->row[row][1],
4967 tab->mat->row[row][0]);
4971 return v;
4974 /* Check if the sequence of variables starting at "pos"
4975 * represents a trivial solution according to "trivial".
4976 * That is, is the result of applying "trivial" to this sequence
4977 * equal to the zero vector?
4979 static isl_bool region_is_trivial(struct isl_tab *tab, int pos,
4980 __isl_keep isl_mat *trivial)
4982 int n, len;
4983 isl_vec *v;
4984 isl_bool is_trivial;
4986 if (!trivial)
4987 return isl_bool_error;
4989 n = isl_mat_rows(trivial);
4990 if (n == 0)
4991 return isl_bool_false;
4993 len = isl_mat_cols(trivial);
4994 v = extract_sample_sequence(tab, pos, len);
4995 v = isl_mat_vec_product(isl_mat_copy(trivial), v);
4996 is_trivial = isl_vec_is_zero(v);
4997 isl_vec_free(v);
4999 return is_trivial;
5002 /* Global internal data for isl_tab_basic_set_non_trivial_lexmin.
5004 * "n_op" is the number of initial coordinates to optimize,
5005 * as passed to isl_tab_basic_set_non_trivial_lexmin.
5006 * "region" is the "n_region"-sized array of regions passed
5007 * to isl_tab_basic_set_non_trivial_lexmin.
5009 * "tab" is the tableau that corresponds to the ILP problem.
5010 * "local" is an array of local data structure, one for each
5011 * (potential) level of the backtracking procedure of
5012 * isl_tab_basic_set_non_trivial_lexmin.
5013 * "v" is a pre-allocated vector that can be used for adding
5014 * constraints to the tableau.
5016 * "sol" contains the best solution found so far.
5017 * It is initialized to a vector of size zero.
5019 struct isl_lexmin_data {
5020 int n_op;
5021 int n_region;
5022 struct isl_trivial_region *region;
5024 struct isl_tab *tab;
5025 struct isl_local_region *local;
5026 isl_vec *v;
5028 isl_vec *sol;
5031 /* Return the index of the first trivial region, "n_region" if all regions
5032 * are non-trivial or -1 in case of error.
5034 static int first_trivial_region(struct isl_lexmin_data *data)
5036 int i;
5038 for (i = 0; i < data->n_region; ++i) {
5039 isl_bool trivial;
5040 trivial = region_is_trivial(data->tab, data->region[i].pos,
5041 data->region[i].trivial);
5042 if (trivial < 0)
5043 return -1;
5044 if (trivial)
5045 return i;
5048 return data->n_region;
5051 /* Check if the solution is optimal, i.e., whether the first
5052 * n_op entries are zero.
5054 static int is_optimal(__isl_keep isl_vec *sol, int n_op)
5056 int i;
5058 for (i = 0; i < n_op; ++i)
5059 if (!isl_int_is_zero(sol->el[1 + i]))
5060 return 0;
5061 return 1;
5064 /* Add constraints to "tab" that ensure that any solution is significantly
5065 * better than that represented by "sol". That is, find the first
5066 * relevant (within first n_op) non-zero coefficient and force it (along
5067 * with all previous coefficients) to be zero.
5068 * If the solution is already optimal (all relevant coefficients are zero),
5069 * then just mark the table as empty.
5070 * "n_zero" is the number of coefficients that have been forced zero
5071 * by previous calls to this function at the same level.
5072 * Return the updated number of forced zero coefficients or -1 on error.
5074 * This function assumes that at least 2 * (n_op - n_zero) more rows and
5075 * at least 2 * (n_op - n_zero) more elements in the constraint array
5076 * are available in the tableau.
5078 static int force_better_solution(struct isl_tab *tab,
5079 __isl_keep isl_vec *sol, int n_op, int n_zero)
5081 int i, n;
5082 isl_ctx *ctx;
5083 isl_vec *v = NULL;
5085 if (!sol)
5086 return -1;
5088 for (i = n_zero; i < n_op; ++i)
5089 if (!isl_int_is_zero(sol->el[1 + i]))
5090 break;
5092 if (i == n_op) {
5093 if (isl_tab_mark_empty(tab) < 0)
5094 return -1;
5095 return n_op;
5098 ctx = isl_vec_get_ctx(sol);
5099 v = isl_vec_alloc(ctx, 1 + tab->n_var);
5100 if (!v)
5101 return -1;
5103 n = i + 1;
5104 for (; i >= n_zero; --i) {
5105 v = isl_vec_clr(v);
5106 isl_int_set_si(v->el[1 + i], -1);
5107 if (add_lexmin_eq(tab, v->el) < 0)
5108 goto error;
5111 isl_vec_free(v);
5112 return n;
5113 error:
5114 isl_vec_free(v);
5115 return -1;
5118 /* Fix triviality direction "dir" of the given region to zero.
5120 * This function assumes that at least two more rows and at least
5121 * two more elements in the constraint array are available in the tableau.
5123 static isl_stat fix_zero(struct isl_tab *tab, struct isl_trivial_region *region,
5124 int dir, struct isl_lexmin_data *data)
5126 int len;
5128 data->v = isl_vec_clr(data->v);
5129 if (!data->v)
5130 return isl_stat_error;
5131 len = isl_mat_cols(region->trivial);
5132 isl_seq_cpy(data->v->el + 1 + region->pos, region->trivial->row[dir],
5133 len);
5134 if (add_lexmin_eq(tab, data->v->el) < 0)
5135 return isl_stat_error;
5137 return isl_stat_ok;
5140 /* This function selects case "side" for non-triviality region "region",
5141 * assuming all the equality constraints have been imposed already.
5142 * In particular, the triviality direction side/2 is made positive
5143 * if side is even and made negative if side is odd.
5145 * This function assumes that at least one more row and at least
5146 * one more element in the constraint array are available in the tableau.
5148 static struct isl_tab *pos_neg(struct isl_tab *tab,
5149 struct isl_trivial_region *region,
5150 int side, struct isl_lexmin_data *data)
5152 int len;
5154 data->v = isl_vec_clr(data->v);
5155 if (!data->v)
5156 goto error;
5157 isl_int_set_si(data->v->el[0], -1);
5158 len = isl_mat_cols(region->trivial);
5159 if (side % 2 == 0)
5160 isl_seq_cpy(data->v->el + 1 + region->pos,
5161 region->trivial->row[side / 2], len);
5162 else
5163 isl_seq_neg(data->v->el + 1 + region->pos,
5164 region->trivial->row[side / 2], len);
5165 return add_lexmin_ineq(tab, data->v->el);
5166 error:
5167 isl_tab_free(tab);
5168 return NULL;
5171 /* Local data at each level of the backtracking procedure of
5172 * isl_tab_basic_set_non_trivial_lexmin.
5174 * "update" is set if a solution has been found in the current case
5175 * of this level, such that a better solution needs to be enforced
5176 * in the next case.
5177 * "n_zero" is the number of initial coordinates that have already
5178 * been forced to be zero at this level.
5179 * "region" is the non-triviality region considered at this level.
5180 * "side" is the index of the current case at this level.
5181 * "n" is the number of triviality directions.
5182 * "snap" is a snapshot of the tableau holding a state that needs
5183 * to be satisfied by all subsequent cases.
5185 struct isl_local_region {
5186 int update;
5187 int n_zero;
5188 int region;
5189 int side;
5190 int n;
5191 struct isl_tab_undo *snap;
5194 /* Initialize the global data structure "data" used while solving
5195 * the ILP problem "bset".
5197 static isl_stat init_lexmin_data(struct isl_lexmin_data *data,
5198 __isl_keep isl_basic_set *bset)
5200 isl_ctx *ctx;
5202 ctx = isl_basic_set_get_ctx(bset);
5204 data->tab = tab_for_lexmin(bset, NULL, 0, 0);
5205 if (!data->tab)
5206 return isl_stat_error;
5208 data->v = isl_vec_alloc(ctx, 1 + data->tab->n_var);
5209 if (!data->v)
5210 return isl_stat_error;
5211 data->local = isl_calloc_array(ctx, struct isl_local_region,
5212 data->n_region);
5213 if (data->n_region && !data->local)
5214 return isl_stat_error;
5216 data->sol = isl_vec_alloc(ctx, 0);
5218 return isl_stat_ok;
5221 /* Mark all outer levels as requiring a better solution
5222 * in the next cases.
5224 static void update_outer_levels(struct isl_lexmin_data *data, int level)
5226 int i;
5228 for (i = 0; i < level; ++i)
5229 data->local[i].update = 1;
5232 /* Initialize "local" to refer to region "region" and
5233 * to initiate processing at this level.
5235 static void init_local_region(struct isl_local_region *local, int region,
5236 struct isl_lexmin_data *data)
5238 local->n = isl_mat_rows(data->region[region].trivial);
5239 local->region = region;
5240 local->side = 0;
5241 local->update = 0;
5242 local->n_zero = 0;
5245 /* What to do next after entering a level of the backtracking procedure.
5247 * error: some error has occurred; abort
5248 * done: an optimal solution has been found; stop search
5249 * backtrack: backtrack to the previous level
5250 * handle: add the constraints for the current level and
5251 * move to the next level
5253 enum isl_next {
5254 isl_next_error = -1,
5255 isl_next_done,
5256 isl_next_backtrack,
5257 isl_next_handle,
5260 /* Have all cases of the current region been considered?
5261 * If there are n directions, then there are 2n cases.
5263 * The constraints in the current tableau are imposed
5264 * in all subsequent cases. This means that if the current
5265 * tableau is empty, then none of those cases should be considered
5266 * anymore and all cases have effectively been considered.
5268 static int finished_all_cases(struct isl_local_region *local,
5269 struct isl_lexmin_data *data)
5271 if (data->tab->empty)
5272 return 1;
5273 return local->side >= 2 * local->n;
5276 /* Enter level "level" of the backtracking search and figure out
5277 * what to do next. "init" is set if the level was entered
5278 * from a higher level and needs to be initialized.
5279 * Otherwise, the level is entered as a result of backtracking and
5280 * the tableau needs to be restored to a position that can
5281 * be used for the next case at this level.
5282 * The snapshot is assumed to have been saved in the previous case,
5283 * before the constraints specific to that case were added.
5285 * In the initialization case, the local region is initialized
5286 * to point to the first violated region.
5287 * If the constraints of all regions are satisfied by the current
5288 * sample of the tableau, then tell the caller to continue looking
5289 * for a better solution or to stop searching if an optimal solution
5290 * has been found.
5292 * If the tableau is empty or if all cases at the current level
5293 * have been considered, then the caller needs to backtrack as well.
5295 static enum isl_next enter_level(int level, int init,
5296 struct isl_lexmin_data *data)
5298 struct isl_local_region *local = &data->local[level];
5300 if (init) {
5301 int r;
5303 data->tab = cut_to_integer_lexmin(data->tab, CUT_ONE);
5304 if (!data->tab)
5305 return isl_next_error;
5306 if (data->tab->empty)
5307 return isl_next_backtrack;
5308 r = first_trivial_region(data);
5309 if (r < 0)
5310 return isl_next_error;
5311 if (r == data->n_region) {
5312 update_outer_levels(data, level);
5313 isl_vec_free(data->sol);
5314 data->sol = isl_tab_get_sample_value(data->tab);
5315 if (!data->sol)
5316 return isl_next_error;
5317 if (is_optimal(data->sol, data->n_op))
5318 return isl_next_done;
5319 return isl_next_backtrack;
5321 if (level >= data->n_region)
5322 isl_die(isl_vec_get_ctx(data->v), isl_error_internal,
5323 "nesting level too deep",
5324 return isl_next_error);
5325 init_local_region(local, r, data);
5326 if (isl_tab_extend_cons(data->tab,
5327 2 * local->n + 2 * data->n_op) < 0)
5328 return isl_next_error;
5329 } else {
5330 if (isl_tab_rollback(data->tab, local->snap) < 0)
5331 return isl_next_error;
5334 if (finished_all_cases(local, data))
5335 return isl_next_backtrack;
5336 return isl_next_handle;
5339 /* If a solution has been found in the previous case at this level
5340 * (marked by local->update being set), then add constraints
5341 * that enforce a better solution in the present and all following cases.
5342 * The constraints only need to be imposed once because they are
5343 * included in the snapshot (taken in pick_side) that will be used in
5344 * subsequent cases.
5346 static isl_stat better_next_side(struct isl_local_region *local,
5347 struct isl_lexmin_data *data)
5349 if (!local->update)
5350 return isl_stat_ok;
5352 local->n_zero = force_better_solution(data->tab,
5353 data->sol, data->n_op, local->n_zero);
5354 if (local->n_zero < 0)
5355 return isl_stat_error;
5357 local->update = 0;
5359 return isl_stat_ok;
5362 /* Add constraints to data->tab that select the current case (local->side)
5363 * at the current level.
5365 * If the linear combinations v should not be zero, then the cases are
5366 * v_0 >= 1
5367 * v_0 <= -1
5368 * v_0 = 0 and v_1 >= 1
5369 * v_0 = 0 and v_1 <= -1
5370 * v_0 = 0 and v_1 = 0 and v_2 >= 1
5371 * v_0 = 0 and v_1 = 0 and v_2 <= -1
5372 * ...
5373 * in this order.
5375 * A snapshot is taken after the equality constraint (if any) has been added
5376 * such that the next case can start off from this position.
5377 * The rollback to this position is performed in enter_level.
5379 static isl_stat pick_side(struct isl_local_region *local,
5380 struct isl_lexmin_data *data)
5382 struct isl_trivial_region *region;
5383 int side, base;
5385 region = &data->region[local->region];
5386 side = local->side;
5387 base = 2 * (side/2);
5389 if (side == base && base >= 2 &&
5390 fix_zero(data->tab, region, base / 2 - 1, data) < 0)
5391 return isl_stat_error;
5393 local->snap = isl_tab_snap(data->tab);
5394 if (isl_tab_push_basis(data->tab) < 0)
5395 return isl_stat_error;
5397 data->tab = pos_neg(data->tab, region, side, data);
5398 if (!data->tab)
5399 return isl_stat_error;
5400 return isl_stat_ok;
5403 /* Free the memory associated to "data".
5405 static void clear_lexmin_data(struct isl_lexmin_data *data)
5407 free(data->local);
5408 isl_vec_free(data->v);
5409 isl_tab_free(data->tab);
5412 /* Return the lexicographically smallest non-trivial solution of the
5413 * given ILP problem.
5415 * All variables are assumed to be non-negative.
5417 * n_op is the number of initial coordinates to optimize.
5418 * That is, once a solution has been found, we will only continue looking
5419 * for solutions that result in significantly better values for those
5420 * initial coordinates. That is, we only continue looking for solutions
5421 * that increase the number of initial zeros in this sequence.
5423 * A solution is non-trivial, if it is non-trivial on each of the
5424 * specified regions. Each region represents a sequence of
5425 * triviality directions on a sequence of variables that starts
5426 * at a given position. A solution is non-trivial on such a region if
5427 * at least one of the triviality directions is non-zero
5428 * on that sequence of variables.
5430 * Whenever a conflict is encountered, all constraints involved are
5431 * reported to the caller through a call to "conflict".
5433 * We perform a simple branch-and-bound backtracking search.
5434 * Each level in the search represents an initially trivial region
5435 * that is forced to be non-trivial.
5436 * At each level we consider 2 * n cases, where n
5437 * is the number of triviality directions.
5438 * In terms of those n directions v_i, we consider the cases
5439 * v_0 >= 1
5440 * v_0 <= -1
5441 * v_0 = 0 and v_1 >= 1
5442 * v_0 = 0 and v_1 <= -1
5443 * v_0 = 0 and v_1 = 0 and v_2 >= 1
5444 * v_0 = 0 and v_1 = 0 and v_2 <= -1
5445 * ...
5446 * in this order.
5448 __isl_give isl_vec *isl_tab_basic_set_non_trivial_lexmin(
5449 __isl_take isl_basic_set *bset, int n_op, int n_region,
5450 struct isl_trivial_region *region,
5451 int (*conflict)(int con, void *user), void *user)
5453 struct isl_lexmin_data data = { n_op, n_region, region };
5454 int level, init;
5456 if (!bset)
5457 return NULL;
5459 if (init_lexmin_data(&data, bset) < 0)
5460 goto error;
5461 data.tab->conflict = conflict;
5462 data.tab->conflict_user = user;
5464 level = 0;
5465 init = 1;
5467 while (level >= 0) {
5468 enum isl_next next;
5469 struct isl_local_region *local = &data.local[level];
5471 next = enter_level(level, init, &data);
5472 if (next < 0)
5473 goto error;
5474 if (next == isl_next_done)
5475 break;
5476 if (next == isl_next_backtrack) {
5477 level--;
5478 init = 0;
5479 continue;
5482 if (better_next_side(local, &data) < 0)
5483 goto error;
5484 if (pick_side(local, &data) < 0)
5485 goto error;
5487 local->side++;
5488 level++;
5489 init = 1;
5492 clear_lexmin_data(&data);
5493 isl_basic_set_free(bset);
5495 return data.sol;
5496 error:
5497 clear_lexmin_data(&data);
5498 isl_basic_set_free(bset);
5499 isl_vec_free(data.sol);
5500 return NULL;
5503 /* Wrapper for a tableau that is used for computing
5504 * the lexicographically smallest rational point of a non-negative set.
5505 * This point is represented by the sample value of "tab",
5506 * unless "tab" is empty.
5508 struct isl_tab_lexmin {
5509 isl_ctx *ctx;
5510 struct isl_tab *tab;
5513 /* Free "tl" and return NULL.
5515 __isl_null isl_tab_lexmin *isl_tab_lexmin_free(__isl_take isl_tab_lexmin *tl)
5517 if (!tl)
5518 return NULL;
5519 isl_ctx_deref(tl->ctx);
5520 isl_tab_free(tl->tab);
5521 free(tl);
5523 return NULL;
5526 /* Construct an isl_tab_lexmin for computing
5527 * the lexicographically smallest rational point in "bset",
5528 * assuming that all variables are non-negative.
5530 __isl_give isl_tab_lexmin *isl_tab_lexmin_from_basic_set(
5531 __isl_take isl_basic_set *bset)
5533 isl_ctx *ctx;
5534 isl_tab_lexmin *tl;
5536 if (!bset)
5537 return NULL;
5539 ctx = isl_basic_set_get_ctx(bset);
5540 tl = isl_calloc_type(ctx, struct isl_tab_lexmin);
5541 if (!tl)
5542 goto error;
5543 tl->ctx = ctx;
5544 isl_ctx_ref(ctx);
5545 tl->tab = tab_for_lexmin(bset, NULL, 0, 0);
5546 isl_basic_set_free(bset);
5547 if (!tl->tab)
5548 return isl_tab_lexmin_free(tl);
5549 return tl;
5550 error:
5551 isl_basic_set_free(bset);
5552 isl_tab_lexmin_free(tl);
5553 return NULL;
5556 /* Return the dimension of the set represented by "tl".
5558 int isl_tab_lexmin_dim(__isl_keep isl_tab_lexmin *tl)
5560 return tl ? tl->tab->n_var : -1;
5563 /* Add the equality with coefficients "eq" to "tl", updating the optimal
5564 * solution if needed.
5565 * The equality is added as two opposite inequality constraints.
5567 __isl_give isl_tab_lexmin *isl_tab_lexmin_add_eq(__isl_take isl_tab_lexmin *tl,
5568 isl_int *eq)
5570 unsigned n_var;
5572 if (!tl || !eq)
5573 return isl_tab_lexmin_free(tl);
5575 if (isl_tab_extend_cons(tl->tab, 2) < 0)
5576 return isl_tab_lexmin_free(tl);
5577 n_var = tl->tab->n_var;
5578 isl_seq_neg(eq, eq, 1 + n_var);
5579 tl->tab = add_lexmin_ineq(tl->tab, eq);
5580 isl_seq_neg(eq, eq, 1 + n_var);
5581 tl->tab = add_lexmin_ineq(tl->tab, eq);
5583 if (!tl->tab)
5584 return isl_tab_lexmin_free(tl);
5586 return tl;
5589 /* Add cuts to "tl" until the sample value reaches an integer value or
5590 * until the result becomes empty.
5592 __isl_give isl_tab_lexmin *isl_tab_lexmin_cut_to_integer(
5593 __isl_take isl_tab_lexmin *tl)
5595 if (!tl)
5596 return NULL;
5597 tl->tab = cut_to_integer_lexmin(tl->tab, CUT_ONE);
5598 if (!tl->tab)
5599 return isl_tab_lexmin_free(tl);
5600 return tl;
5603 /* Return the lexicographically smallest rational point in the basic set
5604 * from which "tl" was constructed.
5605 * If the original input was empty, then return a zero-length vector.
5607 __isl_give isl_vec *isl_tab_lexmin_get_solution(__isl_keep isl_tab_lexmin *tl)
5609 if (!tl)
5610 return NULL;
5611 if (tl->tab->empty)
5612 return isl_vec_alloc(tl->ctx, 0);
5613 else
5614 return isl_tab_get_sample_value(tl->tab);
5617 struct isl_sol_pma {
5618 struct isl_sol sol;
5619 isl_pw_multi_aff *pma;
5620 isl_set *empty;
5623 static void sol_pma_free(struct isl_sol *sol)
5625 struct isl_sol_pma *sol_pma = (struct isl_sol_pma *) sol;
5626 isl_pw_multi_aff_free(sol_pma->pma);
5627 isl_set_free(sol_pma->empty);
5630 /* This function is called for parts of the context where there is
5631 * no solution, with "bset" corresponding to the context tableau.
5632 * Simply add the basic set to the set "empty".
5634 static void sol_pma_add_empty(struct isl_sol_pma *sol,
5635 __isl_take isl_basic_set *bset)
5637 if (!bset || !sol->empty)
5638 goto error;
5640 sol->empty = isl_set_grow(sol->empty, 1);
5641 bset = isl_basic_set_simplify(bset);
5642 bset = isl_basic_set_finalize(bset);
5643 sol->empty = isl_set_add_basic_set(sol->empty, bset);
5644 if (!sol->empty)
5645 sol->sol.error = 1;
5646 return;
5647 error:
5648 isl_basic_set_free(bset);
5649 sol->sol.error = 1;
5652 /* Given a basic set "dom" that represents the context and a tuple of
5653 * affine expressions "maff" defined over this domain, construct
5654 * an isl_pw_multi_aff with a single cell corresponding to "dom" and
5655 * the affine expressions in "maff".
5657 static void sol_pma_add(struct isl_sol_pma *sol,
5658 __isl_take isl_basic_set *dom, __isl_take isl_multi_aff *maff)
5660 isl_pw_multi_aff *pma;
5662 dom = isl_basic_set_simplify(dom);
5663 dom = isl_basic_set_finalize(dom);
5664 pma = isl_pw_multi_aff_alloc(isl_set_from_basic_set(dom), maff);
5665 sol->pma = isl_pw_multi_aff_add_disjoint(sol->pma, pma);
5666 if (!sol->pma)
5667 sol->sol.error = 1;
5670 static void sol_pma_add_empty_wrap(struct isl_sol *sol,
5671 __isl_take isl_basic_set *bset)
5673 sol_pma_add_empty((struct isl_sol_pma *)sol, bset);
5676 static void sol_pma_add_wrap(struct isl_sol *sol,
5677 __isl_take isl_basic_set *dom, __isl_take isl_multi_aff *ma)
5679 sol_pma_add((struct isl_sol_pma *)sol, dom, ma);
5682 /* Construct an isl_sol_pma structure for accumulating the solution.
5683 * If track_empty is set, then we also keep track of the parts
5684 * of the context where there is no solution.
5685 * If max is set, then we are solving a maximization, rather than
5686 * a minimization problem, which means that the variables in the
5687 * tableau have value "M - x" rather than "M + x".
5689 static struct isl_sol *sol_pma_init(__isl_keep isl_basic_map *bmap,
5690 __isl_take isl_basic_set *dom, int track_empty, int max)
5692 struct isl_sol_pma *sol_pma = NULL;
5693 isl_space *space;
5695 if (!bmap)
5696 goto error;
5698 sol_pma = isl_calloc_type(bmap->ctx, struct isl_sol_pma);
5699 if (!sol_pma)
5700 goto error;
5702 sol_pma->sol.free = &sol_pma_free;
5703 if (sol_init(&sol_pma->sol, bmap, dom, max) < 0)
5704 goto error;
5705 sol_pma->sol.add = &sol_pma_add_wrap;
5706 sol_pma->sol.add_empty = track_empty ? &sol_pma_add_empty_wrap : NULL;
5707 space = isl_space_copy(sol_pma->sol.space);
5708 sol_pma->pma = isl_pw_multi_aff_empty(space);
5709 if (!sol_pma->pma)
5710 goto error;
5712 if (track_empty) {
5713 sol_pma->empty = isl_set_alloc_space(isl_basic_set_get_space(dom),
5714 1, ISL_SET_DISJOINT);
5715 if (!sol_pma->empty)
5716 goto error;
5719 isl_basic_set_free(dom);
5720 return &sol_pma->sol;
5721 error:
5722 isl_basic_set_free(dom);
5723 sol_free(&sol_pma->sol);
5724 return NULL;
5727 /* Base case of isl_tab_basic_map_partial_lexopt, after removing
5728 * some obvious symmetries.
5730 * We call basic_map_partial_lexopt_base_sol and extract the results.
5732 static __isl_give isl_pw_multi_aff *basic_map_partial_lexopt_base_pw_multi_aff(
5733 __isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom,
5734 __isl_give isl_set **empty, int max)
5736 isl_pw_multi_aff *result = NULL;
5737 struct isl_sol *sol;
5738 struct isl_sol_pma *sol_pma;
5740 sol = basic_map_partial_lexopt_base_sol(bmap, dom, empty, max,
5741 &sol_pma_init);
5742 if (!sol)
5743 return NULL;
5744 sol_pma = (struct isl_sol_pma *) sol;
5746 result = isl_pw_multi_aff_copy(sol_pma->pma);
5747 if (empty)
5748 *empty = isl_set_copy(sol_pma->empty);
5749 sol_free(&sol_pma->sol);
5750 return result;
5753 /* Given that the last input variable of "maff" represents the minimum
5754 * of some bounds, check whether we need to plug in the expression
5755 * of the minimum.
5757 * In particular, check if the last input variable appears in any
5758 * of the expressions in "maff".
5760 static int need_substitution(__isl_keep isl_multi_aff *maff)
5762 int i;
5763 unsigned pos;
5765 pos = isl_multi_aff_dim(maff, isl_dim_in) - 1;
5767 for (i = 0; i < maff->n; ++i)
5768 if (isl_aff_involves_dims(maff->u.p[i], isl_dim_in, pos, 1))
5769 return 1;
5771 return 0;
5774 /* Given a set of upper bounds on the last "input" variable m,
5775 * construct a piecewise affine expression that selects
5776 * the minimal upper bound to m, i.e.,
5777 * divide the space into cells where one
5778 * of the upper bounds is smaller than all the others and select
5779 * this upper bound on that cell.
5781 * In particular, if there are n bounds b_i, then the result
5782 * consists of n cell, each one of the form
5784 * b_i <= b_j for j > i
5785 * b_i < b_j for j < i
5787 * The affine expression on this cell is
5789 * b_i
5791 static __isl_give isl_pw_aff *set_minimum_pa(__isl_take isl_space *space,
5792 __isl_take isl_mat *var)
5794 int i;
5795 isl_aff *aff = NULL;
5796 isl_basic_set *bset = NULL;
5797 isl_pw_aff *paff = NULL;
5798 isl_space *pw_space;
5799 isl_local_space *ls = NULL;
5801 if (!space || !var)
5802 goto error;
5804 ls = isl_local_space_from_space(isl_space_copy(space));
5805 pw_space = isl_space_copy(space);
5806 pw_space = isl_space_from_domain(pw_space);
5807 pw_space = isl_space_add_dims(pw_space, isl_dim_out, 1);
5808 paff = isl_pw_aff_alloc_size(pw_space, var->n_row);
5810 for (i = 0; i < var->n_row; ++i) {
5811 isl_pw_aff *paff_i;
5813 aff = isl_aff_alloc(isl_local_space_copy(ls));
5814 bset = isl_basic_set_alloc_space(isl_space_copy(space), 0,
5815 0, var->n_row - 1);
5816 if (!aff || !bset)
5817 goto error;
5818 isl_int_set_si(aff->v->el[0], 1);
5819 isl_seq_cpy(aff->v->el + 1, var->row[i], var->n_col);
5820 isl_int_set_si(aff->v->el[1 + var->n_col], 0);
5821 bset = select_minimum(bset, var, i);
5822 paff_i = isl_pw_aff_alloc(isl_set_from_basic_set(bset), aff);
5823 paff = isl_pw_aff_add_disjoint(paff, paff_i);
5826 isl_local_space_free(ls);
5827 isl_space_free(space);
5828 isl_mat_free(var);
5829 return paff;
5830 error:
5831 isl_aff_free(aff);
5832 isl_basic_set_free(bset);
5833 isl_pw_aff_free(paff);
5834 isl_local_space_free(ls);
5835 isl_space_free(space);
5836 isl_mat_free(var);
5837 return NULL;
5840 /* Given a piecewise multi-affine expression of which the last input variable
5841 * is the minimum of the bounds in "cst", plug in the value of the minimum.
5842 * This minimum expression is given in "min_expr_pa".
5843 * The set "min_expr" contains the same information, but in the form of a set.
5844 * The variable is subsequently projected out.
5846 * The implementation is similar to those of "split" and "split_domain".
5847 * If the variable appears in a given expression, then minimum expression
5848 * is plugged in. Otherwise, if the variable appears in the constraints
5849 * and a split is required, then the domain is split. Otherwise, no split
5850 * is performed.
5852 static __isl_give isl_pw_multi_aff *split_domain_pma(
5853 __isl_take isl_pw_multi_aff *opt, __isl_take isl_pw_aff *min_expr_pa,
5854 __isl_take isl_set *min_expr, __isl_take isl_mat *cst)
5856 int n_in;
5857 int i;
5858 isl_space *space;
5859 isl_pw_multi_aff *res;
5861 if (!opt || !min_expr || !cst)
5862 goto error;
5864 n_in = isl_pw_multi_aff_dim(opt, isl_dim_in);
5865 space = isl_pw_multi_aff_get_space(opt);
5866 space = isl_space_drop_dims(space, isl_dim_in, n_in - 1, 1);
5867 res = isl_pw_multi_aff_empty(space);
5869 for (i = 0; i < opt->n; ++i) {
5870 isl_pw_multi_aff *pma;
5872 pma = isl_pw_multi_aff_alloc(isl_set_copy(opt->p[i].set),
5873 isl_multi_aff_copy(opt->p[i].maff));
5874 if (need_substitution(opt->p[i].maff))
5875 pma = isl_pw_multi_aff_substitute(pma,
5876 isl_dim_in, n_in - 1, min_expr_pa);
5877 else {
5878 isl_bool split;
5879 split = need_split_set(opt->p[i].set, cst);
5880 if (split < 0)
5881 pma = isl_pw_multi_aff_free(pma);
5882 else if (split)
5883 pma = isl_pw_multi_aff_intersect_domain(pma,
5884 isl_set_copy(min_expr));
5886 pma = isl_pw_multi_aff_project_out(pma,
5887 isl_dim_in, n_in - 1, 1);
5889 res = isl_pw_multi_aff_add_disjoint(res, pma);
5892 isl_pw_multi_aff_free(opt);
5893 isl_pw_aff_free(min_expr_pa);
5894 isl_set_free(min_expr);
5895 isl_mat_free(cst);
5896 return res;
5897 error:
5898 isl_pw_multi_aff_free(opt);
5899 isl_pw_aff_free(min_expr_pa);
5900 isl_set_free(min_expr);
5901 isl_mat_free(cst);
5902 return NULL;
5905 static __isl_give isl_pw_multi_aff *basic_map_partial_lexopt_pw_multi_aff(
5906 __isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom,
5907 __isl_give isl_set **empty, int max);
5909 /* This function is called from basic_map_partial_lexopt_symm.
5910 * The last variable of "bmap" and "dom" corresponds to the minimum
5911 * of the bounds in "cst". "map_space" is the space of the original
5912 * input relation (of basic_map_partial_lexopt_symm) and "set_space"
5913 * is the space of the original domain.
5915 * We recursively call basic_map_partial_lexopt and then plug in
5916 * the definition of the minimum in the result.
5918 static __isl_give isl_pw_multi_aff *
5919 basic_map_partial_lexopt_symm_core_pw_multi_aff(
5920 __isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom,
5921 __isl_give isl_set **empty, int max, __isl_take isl_mat *cst,
5922 __isl_take isl_space *map_space, __isl_take isl_space *set_space)
5924 isl_pw_multi_aff *opt;
5925 isl_pw_aff *min_expr_pa;
5926 isl_set *min_expr;
5928 min_expr = set_minimum(isl_basic_set_get_space(dom), isl_mat_copy(cst));
5929 min_expr_pa = set_minimum_pa(isl_basic_set_get_space(dom),
5930 isl_mat_copy(cst));
5932 opt = basic_map_partial_lexopt_pw_multi_aff(bmap, dom, empty, max);
5934 if (empty) {
5935 *empty = split(*empty,
5936 isl_set_copy(min_expr), isl_mat_copy(cst));
5937 *empty = isl_set_reset_space(*empty, set_space);
5940 opt = split_domain_pma(opt, min_expr_pa, min_expr, cst);
5941 opt = isl_pw_multi_aff_reset_space(opt, map_space);
5943 return opt;
5946 #undef TYPE
5947 #define TYPE isl_pw_multi_aff
5948 #undef SUFFIX
5949 #define SUFFIX _pw_multi_aff
5950 #include "isl_tab_lexopt_templ.c"