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

source code of polly/lib/External/isl/isl_tab_pip.c