Coverage Report

Created: 2018-12-09 11:54

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