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