1/* 2 * Copyright 2008-2009 Katholieke Universiteit Leuven 3 * Copyright 2013 Ecole Normale Superieure 4 * 5 * Use of this software is governed by the MIT license 6 * 7 * Written by Sven Verdoolaege, K.U.Leuven, Departement 8 * Computerwetenschappen, Celestijnenlaan 200A, B-3001 Leuven, Belgium 9 * and Ecole Normale Superieure, 45 rue d'Ulm, 75230 Paris, France 10 */ 11 12#include <isl_ctx_private.h> 13#include <isl_mat_private.h> 14#include "isl_map_private.h" 15#include "isl_tab.h" 16#include <isl/seq.h> 17#include <isl_config.h> 18 19/* 20 * The implementation of tableaus in this file was inspired by Section 8 21 * of David Detlefs, Greg Nelson and James B. Saxe, "Simplify: a theorem 22 * prover for program checking". 23 */ 24 25struct isl_tab *isl_tab_alloc(struct isl_ctx *ctx, 26 unsigned n_row, unsigned n_var, unsigned M) 27{ 28 int i; 29 struct isl_tab *tab; 30 unsigned off = 2 + M; 31 32 tab = isl_calloc_type(ctx, struct isl_tab); 33 if (!tab) 34 return NULL; 35 tab->mat = isl_mat_alloc(ctx, n_row, off + n_var); 36 if (!tab->mat) 37 goto error; 38 tab->var = isl_alloc_array(ctx, struct isl_tab_var, n_var); 39 if (n_var && !tab->var) 40 goto error; 41 tab->con = isl_alloc_array(ctx, struct isl_tab_var, n_row); 42 if (n_row && !tab->con) 43 goto error; 44 tab->col_var = isl_alloc_array(ctx, int, n_var); 45 if (n_var && !tab->col_var) 46 goto error; 47 tab->row_var = isl_alloc_array(ctx, int, n_row); 48 if (n_row && !tab->row_var) 49 goto error; 50 for (i = 0; i < n_var; ++i) { 51 tab->var[i].index = i; 52 tab->var[i].is_row = 0; 53 tab->var[i].is_nonneg = 0; 54 tab->var[i].is_zero = 0; 55 tab->var[i].is_redundant = 0; 56 tab->var[i].frozen = 0; 57 tab->var[i].negated = 0; 58 tab->col_var[i] = i; 59 } 60 tab->n_row = 0; 61 tab->n_con = 0; 62 tab->n_eq = 0; 63 tab->max_con = n_row; 64 tab->n_col = n_var; 65 tab->n_var = n_var; 66 tab->max_var = n_var; 67 tab->n_param = 0; 68 tab->n_div = 0; 69 tab->n_dead = 0; 70 tab->n_redundant = 0; 71 tab->strict_redundant = 0; 72 tab->need_undo = 0; 73 tab->rational = 0; 74 tab->empty = 0; 75 tab->in_undo = 0; 76 tab->M = M; 77 tab->cone = 0; 78 tab->bottom.type = isl_tab_undo_bottom; 79 tab->bottom.next = NULL; 80 tab->top = &tab->bottom; 81 82 tab->n_zero = 0; 83 tab->n_unbounded = 0; 84 tab->basis = NULL; 85 86 return tab; 87error: 88 isl_tab_free(tab); 89 return NULL; 90} 91 92isl_ctx *isl_tab_get_ctx(struct isl_tab *tab) 93{ 94 return tab ? isl_mat_get_ctx(tab->mat) : NULL; 95} 96 97int isl_tab_extend_cons(struct isl_tab *tab, unsigned n_new) 98{ 99 unsigned off; 100 101 if (!tab) 102 return -1; 103 104 off = 2 + tab->M; 105 106 if (tab->max_con < tab->n_con + n_new) { 107 struct isl_tab_var *con; 108 109 con = isl_realloc_array(tab->mat->ctx, tab->con, 110 struct isl_tab_var, tab->max_con + n_new); 111 if (!con) 112 return -1; 113 tab->con = con; 114 tab->max_con += n_new; 115 } 116 if (tab->mat->n_row < tab->n_row + n_new) { 117 int *row_var; 118 119 tab->mat = isl_mat_extend(tab->mat, 120 tab->n_row + n_new, off + tab->n_col); 121 if (!tab->mat) 122 return -1; 123 row_var = isl_realloc_array(tab->mat->ctx, tab->row_var, 124 int, tab->mat->n_row); 125 if (!row_var) 126 return -1; 127 tab->row_var = row_var; 128 if (tab->row_sign) { 129 enum isl_tab_row_sign *s; 130 s = isl_realloc_array(tab->mat->ctx, tab->row_sign, 131 enum isl_tab_row_sign, tab->mat->n_row); 132 if (!s) 133 return -1; 134 tab->row_sign = s; 135 } 136 } 137 return 0; 138} 139 140/* Make room for at least n_new extra variables. 141 * Return -1 if anything went wrong. 142 */ 143int isl_tab_extend_vars(struct isl_tab *tab, unsigned n_new) 144{ 145 struct isl_tab_var *var; 146 unsigned off = 2 + tab->M; 147 148 if (tab->max_var < tab->n_var + n_new) { 149 var = isl_realloc_array(tab->mat->ctx, tab->var, 150 struct isl_tab_var, tab->n_var + n_new); 151 if (!var) 152 return -1; 153 tab->var = var; 154 tab->max_var += n_new; 155 } 156 157 if (tab->mat->n_col < off + tab->n_col + n_new) { 158 int *p; 159 160 tab->mat = isl_mat_extend(tab->mat, 161 tab->mat->n_row, off + tab->n_col + n_new); 162 if (!tab->mat) 163 return -1; 164 p = isl_realloc_array(tab->mat->ctx, tab->col_var, 165 int, tab->n_col + n_new); 166 if (!p) 167 return -1; 168 tab->col_var = p; 169 } 170 171 return 0; 172} 173 174struct isl_tab *isl_tab_extend(struct isl_tab *tab, unsigned n_new) 175{ 176 if (isl_tab_extend_cons(tab, n_new) >= 0) 177 return tab; 178 179 isl_tab_free(tab); 180 return NULL; 181} 182 183static void free_undo_record(struct isl_tab_undo *undo) 184{ 185 switch (undo->type) { 186 case isl_tab_undo_saved_basis: 187 free(undo->u.col_var); 188 break; 189 default:; 190 } 191 free(undo); 192} 193 194static void free_undo(struct isl_tab *tab) 195{ 196 struct isl_tab_undo *undo, *next; 197 198 for (undo = tab->top; undo && undo != &tab->bottom; undo = next) { 199 next = undo->next; 200 free_undo_record(undo); 201 } 202 tab->top = undo; 203} 204 205void isl_tab_free(struct isl_tab *tab) 206{ 207 if (!tab) 208 return; 209 free_undo(tab); 210 isl_mat_free(tab->mat); 211 isl_vec_free(tab->dual); 212 isl_basic_map_free(tab->bmap); 213 free(tab->var); 214 free(tab->con); 215 free(tab->row_var); 216 free(tab->col_var); 217 free(tab->row_sign); 218 isl_mat_free(tab->samples); 219 free(tab->sample_index); 220 isl_mat_free(tab->basis); 221 free(tab); 222} 223 224struct isl_tab *isl_tab_dup(struct isl_tab *tab) 225{ 226 int i; 227 struct isl_tab *dup; 228 unsigned off; 229 230 if (!tab) 231 return NULL; 232 233 off = 2 + tab->M; 234 dup = isl_calloc_type(tab->mat->ctx, struct isl_tab); 235 if (!dup) 236 return NULL; 237 dup->mat = isl_mat_dup(tab->mat); 238 if (!dup->mat) 239 goto error; 240 dup->var = isl_alloc_array(tab->mat->ctx, struct isl_tab_var, tab->max_var); 241 if (tab->max_var && !dup->var) 242 goto error; 243 for (i = 0; i < tab->n_var; ++i) 244 dup->var[i] = tab->var[i]; 245 dup->con = isl_alloc_array(tab->mat->ctx, struct isl_tab_var, tab->max_con); 246 if (tab->max_con && !dup->con) 247 goto error; 248 for (i = 0; i < tab->n_con; ++i) 249 dup->con[i] = tab->con[i]; 250 dup->col_var = isl_alloc_array(tab->mat->ctx, int, tab->mat->n_col - off); 251 if ((tab->mat->n_col - off) && !dup->col_var) 252 goto error; 253 for (i = 0; i < tab->n_col; ++i) 254 dup->col_var[i] = tab->col_var[i]; 255 dup->row_var = isl_alloc_array(tab->mat->ctx, int, tab->mat->n_row); 256 if (tab->mat->n_row && !dup->row_var) 257 goto error; 258 for (i = 0; i < tab->n_row; ++i) 259 dup->row_var[i] = tab->row_var[i]; 260 if (tab->row_sign) { 261 dup->row_sign = isl_alloc_array(tab->mat->ctx, enum isl_tab_row_sign, 262 tab->mat->n_row); 263 if (tab->mat->n_row && !dup->row_sign) 264 goto error; 265 for (i = 0; i < tab->n_row; ++i) 266 dup->row_sign[i] = tab->row_sign[i]; 267 } 268 if (tab->samples) { 269 dup->samples = isl_mat_dup(tab->samples); 270 if (!dup->samples) 271 goto error; 272 dup->sample_index = isl_alloc_array(tab->mat->ctx, int, 273 tab->samples->n_row); 274 if (tab->samples->n_row && !dup->sample_index) 275 goto error; 276 dup->n_sample = tab->n_sample; 277 dup->n_outside = tab->n_outside; 278 } 279 dup->n_row = tab->n_row; 280 dup->n_con = tab->n_con; 281 dup->n_eq = tab->n_eq; 282 dup->max_con = tab->max_con; 283 dup->n_col = tab->n_col; 284 dup->n_var = tab->n_var; 285 dup->max_var = tab->max_var; 286 dup->n_param = tab->n_param; 287 dup->n_div = tab->n_div; 288 dup->n_dead = tab->n_dead; 289 dup->n_redundant = tab->n_redundant; 290 dup->rational = tab->rational; 291 dup->empty = tab->empty; 292 dup->strict_redundant = 0; 293 dup->need_undo = 0; 294 dup->in_undo = 0; 295 dup->M = tab->M; 296 tab->cone = tab->cone; 297 dup->bottom.type = isl_tab_undo_bottom; 298 dup->bottom.next = NULL; 299 dup->top = &dup->bottom; 300 301 dup->n_zero = tab->n_zero; 302 dup->n_unbounded = tab->n_unbounded; 303 dup->basis = isl_mat_dup(tab->basis); 304 305 return dup; 306error: 307 isl_tab_free(dup); 308 return NULL; 309} 310 311/* Construct the coefficient matrix of the product tableau 312 * of two tableaus. 313 * mat{1,2} is the coefficient matrix of tableau {1,2} 314 * row{1,2} is the number of rows in tableau {1,2} 315 * col{1,2} is the number of columns in tableau {1,2} 316 * off is the offset to the coefficient column (skipping the 317 * denominator, the constant term and the big parameter if any) 318 * r{1,2} is the number of redundant rows in tableau {1,2} 319 * d{1,2} is the number of dead columns in tableau {1,2} 320 * 321 * The order of the rows and columns in the result is as explained 322 * in isl_tab_product. 323 */ 324static struct isl_mat *tab_mat_product(struct isl_mat *mat1, 325 struct isl_mat *mat2, unsigned row1, unsigned row2, 326 unsigned col1, unsigned col2, 327 unsigned off, unsigned r1, unsigned r2, unsigned d1, unsigned d2) 328{ 329 int i; 330 struct isl_mat *prod; 331 unsigned n; 332 333 prod = isl_mat_alloc(mat1->ctx, mat1->n_row + mat2->n_row, 334 off + col1 + col2); 335 if (!prod) 336 return NULL; 337 338 n = 0; 339 for (i = 0; i < r1; ++i) { 340 isl_seq_cpy(prod->row[n + i], mat1->row[i], off + d1); 341 isl_seq_clr(prod->row[n + i] + off + d1, d2); 342 isl_seq_cpy(prod->row[n + i] + off + d1 + d2, 343 mat1->row[i] + off + d1, col1 - d1); 344 isl_seq_clr(prod->row[n + i] + off + col1 + d1, col2 - d2); 345 } 346 347 n += r1; 348 for (i = 0; i < r2; ++i) { 349 isl_seq_cpy(prod->row[n + i], mat2->row[i], off); 350 isl_seq_clr(prod->row[n + i] + off, d1); 351 isl_seq_cpy(prod->row[n + i] + off + d1, 352 mat2->row[i] + off, d2); 353 isl_seq_clr(prod->row[n + i] + off + d1 + d2, col1 - d1); 354 isl_seq_cpy(prod->row[n + i] + off + col1 + d1, 355 mat2->row[i] + off + d2, col2 - d2); 356 } 357 358 n += r2; 359 for (i = 0; i < row1 - r1; ++i) { 360 isl_seq_cpy(prod->row[n + i], mat1->row[r1 + i], off + d1); 361 isl_seq_clr(prod->row[n + i] + off + d1, d2); 362 isl_seq_cpy(prod->row[n + i] + off + d1 + d2, 363 mat1->row[r1 + i] + off + d1, col1 - d1); 364 isl_seq_clr(prod->row[n + i] + off + col1 + d1, col2 - d2); 365 } 366 367 n += row1 - r1; 368 for (i = 0; i < row2 - r2; ++i) { 369 isl_seq_cpy(prod->row[n + i], mat2->row[r2 + i], off); 370 isl_seq_clr(prod->row[n + i] + off, d1); 371 isl_seq_cpy(prod->row[n + i] + off + d1, 372 mat2->row[r2 + i] + off, d2); 373 isl_seq_clr(prod->row[n + i] + off + d1 + d2, col1 - d1); 374 isl_seq_cpy(prod->row[n + i] + off + col1 + d1, 375 mat2->row[r2 + i] + off + d2, col2 - d2); 376 } 377 378 return prod; 379} 380 381/* Update the row or column index of a variable that corresponds 382 * to a variable in the first input tableau. 383 */ 384static void update_index1(struct isl_tab_var *var, 385 unsigned r1, unsigned r2, unsigned d1, unsigned d2) 386{ 387 if (var->index == -1) 388 return; 389 if (var->is_row && var->index >= r1) 390 var->index += r2; 391 if (!var->is_row && var->index >= d1) 392 var->index += d2; 393} 394 395/* Update the row or column index of a variable that corresponds 396 * to a variable in the second input tableau. 397 */ 398static void update_index2(struct isl_tab_var *var, 399 unsigned row1, unsigned col1, 400 unsigned r1, unsigned r2, unsigned d1, unsigned d2) 401{ 402 if (var->index == -1) 403 return; 404 if (var->is_row) { 405 if (var->index < r2) 406 var->index += r1; 407 else 408 var->index += row1; 409 } else { 410 if (var->index < d2) 411 var->index += d1; 412 else 413 var->index += col1; 414 } 415} 416 417/* Create a tableau that represents the Cartesian product of the sets 418 * represented by tableaus tab1 and tab2. 419 * The order of the rows in the product is 420 * - redundant rows of tab1 421 * - redundant rows of tab2 422 * - non-redundant rows of tab1 423 * - non-redundant rows of tab2 424 * The order of the columns is 425 * - denominator 426 * - constant term 427 * - coefficient of big parameter, if any 428 * - dead columns of tab1 429 * - dead columns of tab2 430 * - live columns of tab1 431 * - live columns of tab2 432 * The order of the variables and the constraints is a concatenation 433 * of order in the two input tableaus. 434 */ 435struct isl_tab *isl_tab_product(struct isl_tab *tab1, struct isl_tab *tab2) 436{ 437 int i; 438 struct isl_tab *prod; 439 unsigned off; 440 unsigned r1, r2, d1, d2; 441 442 if (!tab1 || !tab2) 443 return NULL; 444 445 isl_assert(tab1->mat->ctx, tab1->M == tab2->M, return NULL); 446 isl_assert(tab1->mat->ctx, tab1->rational == tab2->rational, return NULL); 447 isl_assert(tab1->mat->ctx, tab1->cone == tab2->cone, return NULL); 448 isl_assert(tab1->mat->ctx, !tab1->row_sign, return NULL); 449 isl_assert(tab1->mat->ctx, !tab2->row_sign, return NULL); 450 isl_assert(tab1->mat->ctx, tab1->n_param == 0, return NULL); 451 isl_assert(tab1->mat->ctx, tab2->n_param == 0, return NULL); 452 isl_assert(tab1->mat->ctx, tab1->n_div == 0, return NULL); 453 isl_assert(tab1->mat->ctx, tab2->n_div == 0, return NULL); 454 455 off = 2 + tab1->M; 456 r1 = tab1->n_redundant; 457 r2 = tab2->n_redundant; 458 d1 = tab1->n_dead; 459 d2 = tab2->n_dead; 460 prod = isl_calloc_type(tab1->mat->ctx, struct isl_tab); 461 if (!prod) 462 return NULL; 463 prod->mat = tab_mat_product(tab1->mat, tab2->mat, 464 tab1->n_row, tab2->n_row, 465 tab1->n_col, tab2->n_col, off, r1, r2, d1, d2); 466 if (!prod->mat) 467 goto error; 468 prod->var = isl_alloc_array(tab1->mat->ctx, struct isl_tab_var, 469 tab1->max_var + tab2->max_var); 470 if ((tab1->max_var + tab2->max_var) && !prod->var) 471 goto error; 472 for (i = 0; i < tab1->n_var; ++i) { 473 prod->var[i] = tab1->var[i]; 474 update_index1(&prod->var[i], r1, r2, d1, d2); 475 } 476 for (i = 0; i < tab2->n_var; ++i) { 477 prod->var[tab1->n_var + i] = tab2->var[i]; 478 update_index2(&prod->var[tab1->n_var + i], 479 tab1->n_row, tab1->n_col, 480 r1, r2, d1, d2); 481 } 482 prod->con = isl_alloc_array(tab1->mat->ctx, struct isl_tab_var, 483 tab1->max_con + tab2->max_con); 484 if ((tab1->max_con + tab2->max_con) && !prod->con) 485 goto error; 486 for (i = 0; i < tab1->n_con; ++i) { 487 prod->con[i] = tab1->con[i]; 488 update_index1(&prod->con[i], r1, r2, d1, d2); 489 } 490 for (i = 0; i < tab2->n_con; ++i) { 491 prod->con[tab1->n_con + i] = tab2->con[i]; 492 update_index2(&prod->con[tab1->n_con + i], 493 tab1->n_row, tab1->n_col, 494 r1, r2, d1, d2); 495 } 496 prod->col_var = isl_alloc_array(tab1->mat->ctx, int, 497 tab1->n_col + tab2->n_col); 498 if ((tab1->n_col + tab2->n_col) && !prod->col_var) 499 goto error; 500 for (i = 0; i < tab1->n_col; ++i) { 501 int pos = i < d1 ? i : i + d2; 502 prod->col_var[pos] = tab1->col_var[i]; 503 } 504 for (i = 0; i < tab2->n_col; ++i) { 505 int pos = i < d2 ? d1 + i : tab1->n_col + i; 506 int t = tab2->col_var[i]; 507 if (t >= 0) 508 t += tab1->n_var; 509 else 510 t -= tab1->n_con; 511 prod->col_var[pos] = t; 512 } 513 prod->row_var = isl_alloc_array(tab1->mat->ctx, int, 514 tab1->mat->n_row + tab2->mat->n_row); 515 if ((tab1->mat->n_row + tab2->mat->n_row) && !prod->row_var) 516 goto error; 517 for (i = 0; i < tab1->n_row; ++i) { 518 int pos = i < r1 ? i : i + r2; 519 prod->row_var[pos] = tab1->row_var[i]; 520 } 521 for (i = 0; i < tab2->n_row; ++i) { 522 int pos = i < r2 ? r1 + i : tab1->n_row + i; 523 int t = tab2->row_var[i]; 524 if (t >= 0) 525 t += tab1->n_var; 526 else 527 t -= tab1->n_con; 528 prod->row_var[pos] = t; 529 } 530 prod->samples = NULL; 531 prod->sample_index = NULL; 532 prod->n_row = tab1->n_row + tab2->n_row; 533 prod->n_con = tab1->n_con + tab2->n_con; 534 prod->n_eq = 0; 535 prod->max_con = tab1->max_con + tab2->max_con; 536 prod->n_col = tab1->n_col + tab2->n_col; 537 prod->n_var = tab1->n_var + tab2->n_var; 538 prod->max_var = tab1->max_var + tab2->max_var; 539 prod->n_param = 0; 540 prod->n_div = 0; 541 prod->n_dead = tab1->n_dead + tab2->n_dead; 542 prod->n_redundant = tab1->n_redundant + tab2->n_redundant; 543 prod->rational = tab1->rational; 544 prod->empty = tab1->empty || tab2->empty; 545 prod->strict_redundant = tab1->strict_redundant || tab2->strict_redundant; 546 prod->need_undo = 0; 547 prod->in_undo = 0; 548 prod->M = tab1->M; 549 prod->cone = tab1->cone; 550 prod->bottom.type = isl_tab_undo_bottom; 551 prod->bottom.next = NULL; 552 prod->top = &prod->bottom; 553 554 prod->n_zero = 0; 555 prod->n_unbounded = 0; 556 prod->basis = NULL; 557 558 return prod; 559error: 560 isl_tab_free(prod); 561 return NULL; 562} 563 564static struct isl_tab_var *var_from_index(struct isl_tab *tab, int i) 565{ 566 if (i >= 0) 567 return &tab->var[i]; 568 else 569 return &tab->con[~i]; 570} 571 572struct isl_tab_var *isl_tab_var_from_row(struct isl_tab *tab, int i) 573{ 574 return var_from_index(tab, tab->row_var[i]); 575} 576 577static struct isl_tab_var *var_from_col(struct isl_tab *tab, int i) 578{ 579 return var_from_index(tab, tab->col_var[i]); 580} 581 582/* Check if there are any upper bounds on column variable "var", 583 * i.e., non-negative rows where var appears with a negative coefficient. 584 * Return 1 if there are no such bounds. 585 */ 586static int max_is_manifestly_unbounded(struct isl_tab *tab, 587 struct isl_tab_var *var) 588{ 589 int i; 590 unsigned off = 2 + tab->M; 591 592 if (var->is_row) 593 return 0; 594 for (i = tab->n_redundant; i < tab->n_row; ++i) { 595 if (!isl_int_is_neg(tab->mat->row[i][off + var->index])) 596 continue; 597 if (isl_tab_var_from_row(tab, i)->is_nonneg) 598 return 0; 599 } 600 return 1; 601} 602 603/* Check if there are any lower bounds on column variable "var", 604 * i.e., non-negative rows where var appears with a positive coefficient. 605 * Return 1 if there are no such bounds. 606 */ 607static int min_is_manifestly_unbounded(struct isl_tab *tab, 608 struct isl_tab_var *var) 609{ 610 int i; 611 unsigned off = 2 + tab->M; 612 613 if (var->is_row) 614 return 0; 615 for (i = tab->n_redundant; i < tab->n_row; ++i) { 616 if (!isl_int_is_pos(tab->mat->row[i][off + var->index])) 617 continue; 618 if (isl_tab_var_from_row(tab, i)->is_nonneg) 619 return 0; 620 } 621 return 1; 622} 623 624static int row_cmp(struct isl_tab *tab, int r1, int r2, int c, isl_int t) 625{ 626 unsigned off = 2 + tab->M; 627 628 if (tab->M) { 629 int s; 630 isl_int_mul(t, tab->mat->row[r1][2], tab->mat->row[r2][off+c]); 631 isl_int_submul(t, tab->mat->row[r2][2], tab->mat->row[r1][off+c]); 632 s = isl_int_sgn(t); 633 if (s) 634 return s; 635 } 636 isl_int_mul(t, tab->mat->row[r1][1], tab->mat->row[r2][off + c]); 637 isl_int_submul(t, tab->mat->row[r2][1], tab->mat->row[r1][off + c]); 638 return isl_int_sgn(t); 639} 640 641/* Given the index of a column "c", return the index of a row 642 * that can be used to pivot the column in, with either an increase 643 * (sgn > 0) or a decrease (sgn < 0) of the corresponding variable. 644 * If "var" is not NULL, then the row returned will be different from 645 * the one associated with "var". 646 * 647 * Each row in the tableau is of the form 648 * 649 * x_r = a_r0 + \sum_i a_ri x_i 650 * 651 * Only rows with x_r >= 0 and with the sign of a_ri opposite to "sgn" 652 * impose any limit on the increase or decrease in the value of x_c 653 * and this bound is equal to a_r0 / |a_rc|. We are therefore looking 654 * for the row with the smallest (most stringent) such bound. 655 * Note that the common denominator of each row drops out of the fraction. 656 * To check if row j has a smaller bound than row r, i.e., 657 * a_j0 / |a_jc| < a_r0 / |a_rc| or a_j0 |a_rc| < a_r0 |a_jc|, 658 * we check if -sign(a_jc) (a_j0 a_rc - a_r0 a_jc) < 0, 659 * where -sign(a_jc) is equal to "sgn". 660 */ 661static int pivot_row(struct isl_tab *tab, 662 struct isl_tab_var *var, int sgn, int c) 663{ 664 int j, r, tsgn; 665 isl_int t; 666 unsigned off = 2 + tab->M; 667 668 isl_int_init(t); 669 r = -1; 670 for (j = tab->n_redundant; j < tab->n_row; ++j) { 671 if (var && j == var->index) 672 continue; 673 if (!isl_tab_var_from_row(tab, j)->is_nonneg) 674 continue; 675 if (sgn * isl_int_sgn(tab->mat->row[j][off + c]) >= 0) 676 continue; 677 if (r < 0) { 678 r = j; 679 continue; 680 } 681 tsgn = sgn * row_cmp(tab, r, j, c, t); 682 if (tsgn < 0 || (tsgn == 0 && 683 tab->row_var[j] < tab->row_var[r])) 684 r = j; 685 } 686 isl_int_clear(t); 687 return r; 688} 689 690/* Find a pivot (row and col) that will increase (sgn > 0) or decrease 691 * (sgn < 0) the value of row variable var. 692 * If not NULL, then skip_var is a row variable that should be ignored 693 * while looking for a pivot row. It is usually equal to var. 694 * 695 * As the given row in the tableau is of the form 696 * 697 * x_r = a_r0 + \sum_i a_ri x_i 698 * 699 * we need to find a column such that the sign of a_ri is equal to "sgn" 700 * (such that an increase in x_i will have the desired effect) or a 701 * column with a variable that may attain negative values. 702 * If a_ri is positive, then we need to move x_i in the same direction 703 * to obtain the desired effect. Otherwise, x_i has to move in the 704 * opposite direction. 705 */ 706static void find_pivot(struct isl_tab *tab, 707 struct isl_tab_var *var, struct isl_tab_var *skip_var, 708 int sgn, int *row, int *col) 709{ 710 int j, r, c; 711 isl_int *tr; 712 713 *row = *col = -1; 714 715 isl_assert(tab->mat->ctx, var->is_row, return); 716 tr = tab->mat->row[var->index] + 2 + tab->M; 717 718 c = -1; 719 for (j = tab->n_dead; j < tab->n_col; ++j) { 720 if (isl_int_is_zero(tr[j])) 721 continue; 722 if (isl_int_sgn(tr[j]) != sgn && 723 var_from_col(tab, j)->is_nonneg) 724 continue; 725 if (c < 0 || tab->col_var[j] < tab->col_var[c]) 726 c = j; 727 } 728 if (c < 0) 729 return; 730 731 sgn *= isl_int_sgn(tr[c]); 732 r = pivot_row(tab, skip_var, sgn, c); 733 *row = r < 0 ? var->index : r; 734 *col = c; 735} 736 737/* Return 1 if row "row" represents an obviously redundant inequality. 738 * This means 739 * - it represents an inequality or a variable 740 * - that is the sum of a non-negative sample value and a positive 741 * combination of zero or more non-negative constraints. 742 */ 743int isl_tab_row_is_redundant(struct isl_tab *tab, int row) 744{ 745 int i; 746 unsigned off = 2 + tab->M; 747 748 if (tab->row_var[row] < 0 && !isl_tab_var_from_row(tab, row)->is_nonneg) 749 return 0; 750 751 if (isl_int_is_neg(tab->mat->row[row][1])) 752 return 0; 753 if (tab->strict_redundant && isl_int_is_zero(tab->mat->row[row][1])) 754 return 0; 755 if (tab->M && isl_int_is_neg(tab->mat->row[row][2])) 756 return 0; 757 758 for (i = tab->n_dead; i < tab->n_col; ++i) { 759 if (isl_int_is_zero(tab->mat->row[row][off + i])) 760 continue; 761 if (tab->col_var[i] >= 0) 762 return 0; 763 if (isl_int_is_neg(tab->mat->row[row][off + i])) 764 return 0; 765 if (!var_from_col(tab, i)->is_nonneg) 766 return 0; 767 } 768 return 1; 769} 770 771static void swap_rows(struct isl_tab *tab, int row1, int row2) 772{ 773 int t; 774 enum isl_tab_row_sign s; 775 776 t = tab->row_var[row1]; 777 tab->row_var[row1] = tab->row_var[row2]; 778 tab->row_var[row2] = t; 779 isl_tab_var_from_row(tab, row1)->index = row1; 780 isl_tab_var_from_row(tab, row2)->index = row2; 781 tab->mat = isl_mat_swap_rows(tab->mat, row1, row2); 782 783 if (!tab->row_sign) 784 return; 785 s = tab->row_sign[row1]; 786 tab->row_sign[row1] = tab->row_sign[row2]; 787 tab->row_sign[row2] = s; 788} 789 790static int push_union(struct isl_tab *tab, 791 enum isl_tab_undo_type type, union isl_tab_undo_val u) WARN_UNUSED; 792static int push_union(struct isl_tab *tab, 793 enum isl_tab_undo_type type, union isl_tab_undo_val u) 794{ 795 struct isl_tab_undo *undo; 796 797 if (!tab) 798 return -1; 799 if (!tab->need_undo) 800 return 0; 801 802 undo = isl_alloc_type(tab->mat->ctx, struct isl_tab_undo); 803 if (!undo) 804 return -1; 805 undo->type = type; 806 undo->u = u; 807 undo->next = tab->top; 808 tab->top = undo; 809 810 return 0; 811} 812 813int isl_tab_push_var(struct isl_tab *tab, 814 enum isl_tab_undo_type type, struct isl_tab_var *var) 815{ 816 union isl_tab_undo_val u; 817 if (var->is_row) 818 u.var_index = tab->row_var[var->index]; 819 else 820 u.var_index = tab->col_var[var->index]; 821 return push_union(tab, type, u); 822} 823 824int isl_tab_push(struct isl_tab *tab, enum isl_tab_undo_type type) 825{ 826 union isl_tab_undo_val u = { 0 }; 827 return push_union(tab, type, u); 828} 829 830/* Push a record on the undo stack describing the current basic 831 * variables, so that the this state can be restored during rollback. 832 */ 833int isl_tab_push_basis(struct isl_tab *tab) 834{ 835 int i; 836 union isl_tab_undo_val u; 837 838 u.col_var = isl_alloc_array(tab->mat->ctx, int, tab->n_col); 839 if (tab->n_col && !u.col_var) 840 return -1; 841 for (i = 0; i < tab->n_col; ++i) 842 u.col_var[i] = tab->col_var[i]; 843 return push_union(tab, isl_tab_undo_saved_basis, u); 844} 845 846int isl_tab_push_callback(struct isl_tab *tab, struct isl_tab_callback *callback) 847{ 848 union isl_tab_undo_val u; 849 u.callback = callback; 850 return push_union(tab, isl_tab_undo_callback, u); 851} 852 853struct isl_tab *isl_tab_init_samples(struct isl_tab *tab) 854{ 855 if (!tab) 856 return NULL; 857 858 tab->n_sample = 0; 859 tab->n_outside = 0; 860 tab->samples = isl_mat_alloc(tab->mat->ctx, 1, 1 + tab->n_var); 861 if (!tab->samples) 862 goto error; 863 tab->sample_index = isl_alloc_array(tab->mat->ctx, int, 1); 864 if (!tab->sample_index) 865 goto error; 866 return tab; 867error: 868 isl_tab_free(tab); 869 return NULL; 870} 871 872struct isl_tab *isl_tab_add_sample(struct isl_tab *tab, 873 __isl_take isl_vec *sample) 874{ 875 if (!tab || !sample) 876 goto error; 877 878 if (tab->n_sample + 1 > tab->samples->n_row) { 879 int *t = isl_realloc_array(tab->mat->ctx, 880 tab->sample_index, int, tab->n_sample + 1); 881 if (!t) 882 goto error; 883 tab->sample_index = t; 884 } 885 886 tab->samples = isl_mat_extend(tab->samples, 887 tab->n_sample + 1, tab->samples->n_col); 888 if (!tab->samples) 889 goto error; 890 891 isl_seq_cpy(tab->samples->row[tab->n_sample], sample->el, sample->size); 892 isl_vec_free(sample); 893 tab->sample_index[tab->n_sample] = tab->n_sample; 894 tab->n_sample++; 895 896 return tab; 897error: 898 isl_vec_free(sample); 899 isl_tab_free(tab); 900 return NULL; 901} 902 903struct isl_tab *isl_tab_drop_sample(struct isl_tab *tab, int s) 904{ 905 if (s != tab->n_outside) { 906 int t = tab->sample_index[tab->n_outside]; 907 tab->sample_index[tab->n_outside] = tab->sample_index[s]; 908 tab->sample_index[s] = t; 909 isl_mat_swap_rows(tab->samples, tab->n_outside, s); 910 } 911 tab->n_outside++; 912 if (isl_tab_push(tab, isl_tab_undo_drop_sample) < 0) { 913 isl_tab_free(tab); 914 return NULL; 915 } 916 917 return tab; 918} 919 920/* Record the current number of samples so that we can remove newer 921 * samples during a rollback. 922 */ 923int isl_tab_save_samples(struct isl_tab *tab) 924{ 925 union isl_tab_undo_val u; 926 927 if (!tab) 928 return -1; 929 930 u.n = tab->n_sample; 931 return push_union(tab, isl_tab_undo_saved_samples, u); 932} 933 934/* Mark row with index "row" as being redundant. 935 * If we may need to undo the operation or if the row represents 936 * a variable of the original problem, the row is kept, 937 * but no longer considered when looking for a pivot row. 938 * Otherwise, the row is simply removed. 939 * 940 * The row may be interchanged with some other row. If it 941 * is interchanged with a later row, return 1. Otherwise return 0. 942 * If the rows are checked in order in the calling function, 943 * then a return value of 1 means that the row with the given 944 * row number may now contain a different row that hasn't been checked yet. 945 */ 946int isl_tab_mark_redundant(struct isl_tab *tab, int row) 947{ 948 struct isl_tab_var *var = isl_tab_var_from_row(tab, row); 949 var->is_redundant = 1; 950 isl_assert(tab->mat->ctx, row >= tab->n_redundant, return -1); 951 if (tab->preserve || tab->need_undo || tab->row_var[row] >= 0) { 952 if (tab->row_var[row] >= 0 && !var->is_nonneg) { 953 var->is_nonneg = 1; 954 if (isl_tab_push_var(tab, isl_tab_undo_nonneg, var) < 0) 955 return -1; 956 } 957 if (row != tab->n_redundant) 958 swap_rows(tab, row, tab->n_redundant); 959 tab->n_redundant++; 960 return isl_tab_push_var(tab, isl_tab_undo_redundant, var); 961 } else { 962 if (row != tab->n_row - 1) 963 swap_rows(tab, row, tab->n_row - 1); 964 isl_tab_var_from_row(tab, tab->n_row - 1)->index = -1; 965 tab->n_row--; 966 return 1; 967 } 968} 969 970int isl_tab_mark_empty(struct isl_tab *tab) 971{ 972 if (!tab) 973 return -1; 974 if (!tab->empty && tab->need_undo) 975 if (isl_tab_push(tab, isl_tab_undo_empty) < 0) 976 return -1; 977 tab->empty = 1; 978 return 0; 979} 980 981int isl_tab_freeze_constraint(struct isl_tab *tab, int con) 982{ 983 struct isl_tab_var *var; 984 985 if (!tab) 986 return -1; 987 988 var = &tab->con[con]; 989 if (var->frozen) 990 return 0; 991 if (var->index < 0) 992 return 0; 993 var->frozen = 1; 994 995 if (tab->need_undo) 996 return isl_tab_push_var(tab, isl_tab_undo_freeze, var); 997 998 return 0; 999} 1000 1001/* Update the rows signs after a pivot of "row" and "col", with "row_sgn" 1002 * the original sign of the pivot element. 1003 * We only keep track of row signs during PILP solving and in this case 1004 * we only pivot a row with negative sign (meaning the value is always 1005 * non-positive) using a positive pivot element. 1006 * 1007 * For each row j, the new value of the parametric constant is equal to 1008 * 1009 * a_j0 - a_jc a_r0/a_rc 1010 * 1011 * where a_j0 is the original parametric constant, a_rc is the pivot element, 1012 * a_r0 is the parametric constant of the pivot row and a_jc is the 1013 * pivot column entry of the row j. 1014 * Since a_r0 is non-positive and a_rc is positive, the sign of row j 1015 * remains the same if a_jc has the same sign as the row j or if 1016 * a_jc is zero. In all other cases, we reset the sign to "unknown". 1017 */ 1018static void update_row_sign(struct isl_tab *tab, int row, int col, int row_sgn) 1019{ 1020 int i; 1021 struct isl_mat *mat = tab->mat; 1022 unsigned off = 2 + tab->M; 1023 1024 if (!tab->row_sign) 1025 return; 1026 1027 if (tab->row_sign[row] == 0) 1028 return; 1029 isl_assert(mat->ctx, row_sgn > 0, return); 1030 isl_assert(mat->ctx, tab->row_sign[row] == isl_tab_row_neg, return); 1031 tab->row_sign[row] = isl_tab_row_pos; 1032 for (i = 0; i < tab->n_row; ++i) { 1033 int s; 1034 if (i == row) 1035 continue; 1036 s = isl_int_sgn(mat->row[i][off + col]); 1037 if (!s) 1038 continue; 1039 if (!tab->row_sign[i]) 1040 continue; 1041 if (s < 0 && tab->row_sign[i] == isl_tab_row_neg) 1042 continue; 1043 if (s > 0 && tab->row_sign[i] == isl_tab_row_pos) 1044 continue; 1045 tab->row_sign[i] = isl_tab_row_unknown; 1046 } 1047} 1048 1049/* Given a row number "row" and a column number "col", pivot the tableau 1050 * such that the associated variables are interchanged. 1051 * The given row in the tableau expresses 1052 * 1053 * x_r = a_r0 + \sum_i a_ri x_i 1054 * 1055 * or 1056 * 1057 * x_c = 1/a_rc x_r - a_r0/a_rc + sum_{i \ne r} -a_ri/a_rc 1058 * 1059 * Substituting this equality into the other rows 1060 * 1061 * x_j = a_j0 + \sum_i a_ji x_i 1062 * 1063 * with a_jc \ne 0, we obtain 1064 * 1065 * x_j = a_jc/a_rc x_r + a_j0 - a_jc a_r0/a_rc + sum a_ji - a_jc a_ri/a_rc 1066 * 1067 * The tableau 1068 * 1069 * n_rc/d_r n_ri/d_r 1070 * n_jc/d_j n_ji/d_j 1071 * 1072 * where i is any other column and j is any other row, 1073 * is therefore transformed into 1074 * 1075 * s(n_rc)d_r/|n_rc| -s(n_rc)n_ri/|n_rc| 1076 * s(n_rc)d_r n_jc/(|n_rc| d_j) (n_ji |n_rc| - s(n_rc)n_jc n_ri)/(|n_rc| d_j) 1077 * 1078 * The transformation is performed along the following steps 1079 * 1080 * d_r/n_rc n_ri/n_rc 1081 * n_jc/d_j n_ji/d_j 1082 * 1083 * s(n_rc)d_r/|n_rc| -s(n_rc)n_ri/|n_rc| 1084 * n_jc/d_j n_ji/d_j 1085 * 1086 * s(n_rc)d_r/|n_rc| -s(n_rc)n_ri/|n_rc| 1087 * n_jc/(|n_rc| d_j) n_ji/(|n_rc| d_j) 1088 * 1089 * s(n_rc)d_r/|n_rc| -s(n_rc)n_ri/|n_rc| 1090 * n_jc/(|n_rc| d_j) (n_ji |n_rc|)/(|n_rc| d_j) 1091 * 1092 * s(n_rc)d_r/|n_rc| -s(n_rc)n_ri/|n_rc| 1093 * n_jc/(|n_rc| d_j) (n_ji |n_rc| - s(n_rc)n_jc n_ri)/(|n_rc| d_j) 1094 * 1095 * s(n_rc)d_r/|n_rc| -s(n_rc)n_ri/|n_rc| 1096 * s(n_rc)d_r n_jc/(|n_rc| d_j) (n_ji |n_rc| - s(n_rc)n_jc n_ri)/(|n_rc| d_j) 1097 * 1098 */ 1099int isl_tab_pivot(struct isl_tab *tab, int row, int col) 1100{ 1101 int i, j; 1102 int sgn; 1103 int t; 1104 struct isl_mat *mat = tab->mat; 1105 struct isl_tab_var *var; 1106 unsigned off = 2 + tab->M; 1107 1108 if (tab->mat->ctx->abort) { 1109 isl_ctx_set_error(tab->mat->ctx, isl_error_abort); 1110 return -1; 1111 } 1112 1113 isl_int_swap(mat->row[row][0], mat->row[row][off + col]); 1114 sgn = isl_int_sgn(mat->row[row][0]); 1115 if (sgn < 0) { 1116 isl_int_neg(mat->row[row][0], mat->row[row][0]); 1117 isl_int_neg(mat->row[row][off + col], mat->row[row][off + col]); 1118 } else 1119 for (j = 0; j < off - 1 + tab->n_col; ++j) { 1120 if (j == off - 1 + col) 1121 continue; 1122 isl_int_neg(mat->row[row][1 + j], mat->row[row][1 + j]); 1123 } 1124 if (!isl_int_is_one(mat->row[row][0])) 1125 isl_seq_normalize(mat->ctx, mat->row[row], off + tab->n_col); 1126 for (i = 0; i < tab->n_row; ++i) { 1127 if (i == row) 1128 continue; 1129 if (isl_int_is_zero(mat->row[i][off + col])) 1130 continue; 1131 isl_int_mul(mat->row[i][0], mat->row[i][0], mat->row[row][0]); 1132 for (j = 0; j < off - 1 + tab->n_col; ++j) { 1133 if (j == off - 1 + col) 1134 continue; 1135 isl_int_mul(mat->row[i][1 + j], 1136 mat->row[i][1 + j], mat->row[row][0]); 1137 isl_int_addmul(mat->row[i][1 + j], 1138 mat->row[i][off + col], mat->row[row][1 + j]); 1139 } 1140 isl_int_mul(mat->row[i][off + col], 1141 mat->row[i][off + col], mat->row[row][off + col]); 1142 if (!isl_int_is_one(mat->row[i][0])) 1143 isl_seq_normalize(mat->ctx, mat->row[i], off + tab->n_col); 1144 } 1145 t = tab->row_var[row]; 1146 tab->row_var[row] = tab->col_var[col]; 1147 tab->col_var[col] = t; 1148 var = isl_tab_var_from_row(tab, row); 1149 var->is_row = 1; 1150 var->index = row; 1151 var = var_from_col(tab, col); 1152 var->is_row = 0; 1153 var->index = col; 1154 update_row_sign(tab, row, col, sgn); 1155 if (tab->in_undo) 1156 return 0; 1157 for (i = tab->n_redundant; i < tab->n_row; ++i) { 1158 if (isl_int_is_zero(mat->row[i][off + col])) 1159 continue; 1160 if (!isl_tab_var_from_row(tab, i)->frozen && 1161 isl_tab_row_is_redundant(tab, i)) { 1162 int redo = isl_tab_mark_redundant(tab, i); 1163 if (redo < 0) 1164 return -1; 1165 if (redo) 1166 --i; 1167 } 1168 } 1169 return 0; 1170} 1171 1172/* If "var" represents a column variable, then pivot is up (sgn > 0) 1173 * or down (sgn < 0) to a row. The variable is assumed not to be 1174 * unbounded in the specified direction. 1175 * If sgn = 0, then the variable is unbounded in both directions, 1176 * and we pivot with any row we can find. 1177 */ 1178static int to_row(struct isl_tab *tab, struct isl_tab_var *var, int sign) WARN_UNUSED; 1179static int to_row(struct isl_tab *tab, struct isl_tab_var *var, int sign) 1180{ 1181 int r; 1182 unsigned off = 2 + tab->M; 1183 1184 if (var->is_row) 1185 return 0; 1186 1187 if (sign == 0) { 1188 for (r = tab->n_redundant; r < tab->n_row; ++r) 1189 if (!isl_int_is_zero(tab->mat->row[r][off+var->index])) 1190 break; 1191 isl_assert(tab->mat->ctx, r < tab->n_row, return -1); 1192 } else { 1193 r = pivot_row(tab, NULL, sign, var->index); 1194 isl_assert(tab->mat->ctx, r >= 0, return -1); 1195 } 1196 1197 return isl_tab_pivot(tab, r, var->index); 1198} 1199 1200/* Check whether all variables that are marked as non-negative 1201 * also have a non-negative sample value. This function is not 1202 * called from the current code but is useful during debugging. 1203 */ 1204static void check_table(struct isl_tab *tab) __attribute__ ((unused)); 1205static void check_table(struct isl_tab *tab) 1206{ 1207 int i; 1208 1209 if (tab->empty) 1210 return; 1211 for (i = tab->n_redundant; i < tab->n_row; ++i) { 1212 struct isl_tab_var *var; 1213 var = isl_tab_var_from_row(tab, i); 1214 if (!var->is_nonneg) 1215 continue; 1216 if (tab->M) { 1217 isl_assert(tab->mat->ctx, 1218 !isl_int_is_neg(tab->mat->row[i][2]), abort()); 1219 if (isl_int_is_pos(tab->mat->row[i][2])) 1220 continue; 1221 } 1222 isl_assert(tab->mat->ctx, !isl_int_is_neg(tab->mat->row[i][1]), 1223 abort()); 1224 } 1225} 1226 1227/* Return the sign of the maximal value of "var". 1228 * If the sign is not negative, then on return from this function, 1229 * the sample value will also be non-negative. 1230 * 1231 * If "var" is manifestly unbounded wrt positive values, we are done. 1232 * Otherwise, we pivot the variable up to a row if needed 1233 * Then we continue pivoting down until either 1234 * - no more down pivots can be performed 1235 * - the sample value is positive 1236 * - the variable is pivoted into a manifestly unbounded column 1237 */ 1238static int sign_of_max(struct isl_tab *tab, struct isl_tab_var *var) 1239{ 1240 int row, col; 1241 1242 if (max_is_manifestly_unbounded(tab, var)) 1243 return 1; 1244 if (to_row(tab, var, 1) < 0) 1245 return -2; 1246 while (!isl_int_is_pos(tab->mat->row[var->index][1])) { 1247 find_pivot(tab, var, var, 1, &row, &col); 1248 if (row == -1) 1249 return isl_int_sgn(tab->mat->row[var->index][1]); 1250 if (isl_tab_pivot(tab, row, col) < 0) 1251 return -2; 1252 if (!var->is_row) /* manifestly unbounded */ 1253 return 1; 1254 } 1255 return 1; 1256} 1257 1258int isl_tab_sign_of_max(struct isl_tab *tab, int con) 1259{ 1260 struct isl_tab_var *var; 1261 1262 if (!tab) 1263 return -2; 1264 1265 var = &tab->con[con]; 1266 isl_assert(tab->mat->ctx, !var->is_redundant, return -2); 1267 isl_assert(tab->mat->ctx, !var->is_zero, return -2); 1268 1269 return sign_of_max(tab, var); 1270} 1271 1272static int row_is_neg(struct isl_tab *tab, int row) 1273{ 1274 if (!tab->M) 1275 return isl_int_is_neg(tab->mat->row[row][1]); 1276 if (isl_int_is_pos(tab->mat->row[row][2])) 1277 return 0; 1278 if (isl_int_is_neg(tab->mat->row[row][2])) 1279 return 1; 1280 return isl_int_is_neg(tab->mat->row[row][1]); 1281} 1282 1283static int row_sgn(struct isl_tab *tab, int row) 1284{ 1285 if (!tab->M) 1286 return isl_int_sgn(tab->mat->row[row][1]); 1287 if (!isl_int_is_zero(tab->mat->row[row][2])) 1288 return isl_int_sgn(tab->mat->row[row][2]); 1289 else 1290 return isl_int_sgn(tab->mat->row[row][1]); 1291} 1292 1293/* Perform pivots until the row variable "var" has a non-negative 1294 * sample value or until no more upward pivots can be performed. 1295 * Return the sign of the sample value after the pivots have been 1296 * performed. 1297 */ 1298static int restore_row(struct isl_tab *tab, struct isl_tab_var *var) 1299{ 1300 int row, col; 1301 1302 while (row_is_neg(tab, var->index)) { 1303 find_pivot(tab, var, var, 1, &row, &col); 1304 if (row == -1) 1305 break; 1306 if (isl_tab_pivot(tab, row, col) < 0) 1307 return -2; 1308 if (!var->is_row) /* manifestly unbounded */ 1309 return 1; 1310 } 1311 return row_sgn(tab, var->index); 1312} 1313 1314/* Perform pivots until we are sure that the row variable "var" 1315 * can attain non-negative values. After return from this 1316 * function, "var" is still a row variable, but its sample 1317 * value may not be non-negative, even if the function returns 1. 1318 */ 1319static int at_least_zero(struct isl_tab *tab, struct isl_tab_var *var) 1320{ 1321 int row, col; 1322 1323 while (isl_int_is_neg(tab->mat->row[var->index][1])) { 1324 find_pivot(tab, var, var, 1, &row, &col); 1325 if (row == -1) 1326 break; 1327 if (row == var->index) /* manifestly unbounded */ 1328 return 1; 1329 if (isl_tab_pivot(tab, row, col) < 0) 1330 return -1; 1331 } 1332 return !isl_int_is_neg(tab->mat->row[var->index][1]); 1333} 1334 1335/* Return a negative value if "var" can attain negative values. 1336 * Return a non-negative value otherwise. 1337 * 1338 * If "var" is manifestly unbounded wrt negative values, we are done. 1339 * Otherwise, if var is in a column, we can pivot it down to a row. 1340 * Then we continue pivoting down until either 1341 * - the pivot would result in a manifestly unbounded column 1342 * => we don't perform the pivot, but simply return -1 1343 * - no more down pivots can be performed 1344 * - the sample value is negative 1345 * If the sample value becomes negative and the variable is supposed 1346 * to be nonnegative, then we undo the last pivot. 1347 * However, if the last pivot has made the pivoting variable 1348 * obviously redundant, then it may have moved to another row. 1349 * In that case we look for upward pivots until we reach a non-negative 1350 * value again. 1351 */ 1352static int sign_of_min(struct isl_tab *tab, struct isl_tab_var *var) 1353{ 1354 int row, col; 1355 struct isl_tab_var *pivot_var = NULL; 1356 1357 if (min_is_manifestly_unbounded(tab, var)) 1358 return -1; 1359 if (!var->is_row) { 1360 col = var->index; 1361 row = pivot_row(tab, NULL, -1, col); 1362 pivot_var = var_from_col(tab, col); 1363 if (isl_tab_pivot(tab, row, col) < 0) 1364 return -2; 1365 if (var->is_redundant) 1366 return 0; 1367 if (isl_int_is_neg(tab->mat->row[var->index][1])) { 1368 if (var->is_nonneg) { 1369 if (!pivot_var->is_redundant && 1370 pivot_var->index == row) { 1371 if (isl_tab_pivot(tab, row, col) < 0) 1372 return -2; 1373 } else 1374 if (restore_row(tab, var) < -1) 1375 return -2; 1376 } 1377 return -1; 1378 } 1379 } 1380 if (var->is_redundant) 1381 return 0; 1382 while (!isl_int_is_neg(tab->mat->row[var->index][1])) { 1383 find_pivot(tab, var, var, -1, &row, &col); 1384 if (row == var->index) 1385 return -1; 1386 if (row == -1) 1387 return isl_int_sgn(tab->mat->row[var->index][1]); 1388 pivot_var = var_from_col(tab, col); 1389 if (isl_tab_pivot(tab, row, col) < 0) 1390 return -2; 1391 if (var->is_redundant) 1392 return 0; 1393 } 1394 if (pivot_var && var->is_nonneg) { 1395 /* pivot back to non-negative value */ 1396 if (!pivot_var->is_redundant && pivot_var->index == row) { 1397 if (isl_tab_pivot(tab, row, col) < 0) 1398 return -2; 1399 } else 1400 if (restore_row(tab, var) < -1) 1401 return -2; 1402 } 1403 return -1; 1404} 1405 1406static int row_at_most_neg_one(struct isl_tab *tab, int row) 1407{ 1408 if (tab->M) { 1409 if (isl_int_is_pos(tab->mat->row[row][2])) 1410 return 0; 1411 if (isl_int_is_neg(tab->mat->row[row][2])) 1412 return 1; 1413 } 1414 return isl_int_is_neg(tab->mat->row[row][1]) && 1415 isl_int_abs_ge(tab->mat->row[row][1], 1416 tab->mat->row[row][0]); 1417} 1418 1419/* Return 1 if "var" can attain values <= -1. 1420 * Return 0 otherwise. 1421 * 1422 * The sample value of "var" is assumed to be non-negative when the 1423 * the function is called. If 1 is returned then the constraint 1424 * is not redundant and the sample value is made non-negative again before 1425 * the function returns. 1426 */ 1427int isl_tab_min_at_most_neg_one(struct isl_tab *tab, struct isl_tab_var *var) 1428{ 1429 int row, col; 1430 struct isl_tab_var *pivot_var; 1431 1432 if (min_is_manifestly_unbounded(tab, var)) 1433 return 1; 1434 if (!var->is_row) { 1435 col = var->index; 1436 row = pivot_row(tab, NULL, -1, col); 1437 pivot_var = var_from_col(tab, col); 1438 if (isl_tab_pivot(tab, row, col) < 0) 1439 return -1; 1440 if (var->is_redundant) 1441 return 0; 1442 if (row_at_most_neg_one(tab, var->index)) { 1443 if (var->is_nonneg) { 1444 if (!pivot_var->is_redundant && 1445 pivot_var->index == row) { 1446 if (isl_tab_pivot(tab, row, col) < 0) 1447 return -1; 1448 } else 1449 if (restore_row(tab, var) < -1) 1450 return -1; 1451 } 1452 return 1; 1453 } 1454 } 1455 if (var->is_redundant) 1456 return 0; 1457 do { 1458 find_pivot(tab, var, var, -1, &row, &col); 1459 if (row == var->index) { 1460 if (restore_row(tab, var) < -1) 1461 return -1; 1462 return 1; 1463 } 1464 if (row == -1) 1465 return 0; 1466 pivot_var = var_from_col(tab, col); 1467 if (isl_tab_pivot(tab, row, col) < 0) 1468 return -1; 1469 if (var->is_redundant) 1470 return 0; 1471 } while (!row_at_most_neg_one(tab, var->index)); 1472 if (var->is_nonneg) { 1473 /* pivot back to non-negative value */ 1474 if (!pivot_var->is_redundant && pivot_var->index == row) 1475 if (isl_tab_pivot(tab, row, col) < 0) 1476 return -1; 1477 if (restore_row(tab, var) < -1) 1478 return -1; 1479 } 1480 return 1; 1481} 1482 1483/* Return 1 if "var" can attain values >= 1. 1484 * Return 0 otherwise. 1485 */ 1486static int at_least_one(struct isl_tab *tab, struct isl_tab_var *var) 1487{ 1488 int row, col; 1489 isl_int *r; 1490 1491 if (max_is_manifestly_unbounded(tab, var)) 1492 return 1; 1493 if (to_row(tab, var, 1) < 0) 1494 return -1; 1495 r = tab->mat->row[var->index]; 1496 while (isl_int_lt(r[1], r[0])) { 1497 find_pivot(tab, var, var, 1, &row, &col); 1498 if (row == -1) 1499 return isl_int_ge(r[1], r[0]); 1500 if (row == var->index) /* manifestly unbounded */ 1501 return 1; 1502 if (isl_tab_pivot(tab, row, col) < 0) 1503 return -1; 1504 } 1505 return 1; 1506} 1507 1508static void swap_cols(struct isl_tab *tab, int col1, int col2) 1509{ 1510 int t; 1511 unsigned off = 2 + tab->M; 1512 t = tab->col_var[col1]; 1513 tab->col_var[col1] = tab->col_var[col2]; 1514 tab->col_var[col2] = t; 1515 var_from_col(tab, col1)->index = col1; 1516 var_from_col(tab, col2)->index = col2; 1517 tab->mat = isl_mat_swap_cols(tab->mat, off + col1, off + col2); 1518} 1519 1520/* Mark column with index "col" as representing a zero variable. 1521 * If we may need to undo the operation the column is kept, 1522 * but no longer considered. 1523 * Otherwise, the column is simply removed. 1524 * 1525 * The column may be interchanged with some other column. If it 1526 * is interchanged with a later column, return 1. Otherwise return 0. 1527 * If the columns are checked in order in the calling function, 1528 * then a return value of 1 means that the column with the given 1529 * column number may now contain a different column that 1530 * hasn't been checked yet. 1531 */ 1532int isl_tab_kill_col(struct isl_tab *tab, int col) 1533{ 1534 var_from_col(tab, col)->is_zero = 1; 1535 if (tab->need_undo) { 1536 if (isl_tab_push_var(tab, isl_tab_undo_zero, 1537 var_from_col(tab, col)) < 0) 1538 return -1; 1539 if (col != tab->n_dead) 1540 swap_cols(tab, col, tab->n_dead); 1541 tab->n_dead++; 1542 return 0; 1543 } else { 1544 if (col != tab->n_col - 1) 1545 swap_cols(tab, col, tab->n_col - 1); 1546 var_from_col(tab, tab->n_col - 1)->index = -1; 1547 tab->n_col--; 1548 return 1; 1549 } 1550} 1551 1552static int row_is_manifestly_non_integral(struct isl_tab *tab, int row) 1553{ 1554 unsigned off = 2 + tab->M; 1555 1556 if (tab->M && !isl_int_eq(tab->mat->row[row][2], 1557 tab->mat->row[row][0])) 1558 return 0; 1559 if (isl_seq_first_non_zero(tab->mat->row[row] + off + tab->n_dead, 1560 tab->n_col - tab->n_dead) != -1) 1561 return 0; 1562 1563 return !isl_int_is_divisible_by(tab->mat->row[row][1], 1564 tab->mat->row[row][0]); 1565} 1566 1567/* For integer tableaus, check if any of the coordinates are stuck 1568 * at a non-integral value. 1569 */ 1570static int tab_is_manifestly_empty(struct isl_tab *tab) 1571{ 1572 int i; 1573 1574 if (tab->empty) 1575 return 1; 1576 if (tab->rational) 1577 return 0; 1578 1579 for (i = 0; i < tab->n_var; ++i) { 1580 if (!tab->var[i].is_row) 1581 continue; 1582 if (row_is_manifestly_non_integral(tab, tab->var[i].index)) 1583 return 1; 1584 } 1585 1586 return 0; 1587} 1588 1589/* Row variable "var" is non-negative and cannot attain any values 1590 * larger than zero. This means that the coefficients of the unrestricted 1591 * column variables are zero and that the coefficients of the non-negative 1592 * column variables are zero or negative. 1593 * Each of the non-negative variables with a negative coefficient can 1594 * then also be written as the negative sum of non-negative variables 1595 * and must therefore also be zero. 1596 */ 1597static int close_row(struct isl_tab *tab, struct isl_tab_var *var) WARN_UNUSED; 1598static int close_row(struct isl_tab *tab, struct isl_tab_var *var) 1599{ 1600 int j; 1601 struct isl_mat *mat = tab->mat; 1602 unsigned off = 2 + tab->M; 1603 1604 isl_assert(tab->mat->ctx, var->is_nonneg, return -1); 1605 var->is_zero = 1; 1606 if (tab->need_undo) 1607 if (isl_tab_push_var(tab, isl_tab_undo_zero, var) < 0) 1608 return -1; 1609 for (j = tab->n_dead; j < tab->n_col; ++j) { 1610 int recheck; 1611 if (isl_int_is_zero(mat->row[var->index][off + j])) 1612 continue; 1613 isl_assert(tab->mat->ctx, 1614 isl_int_is_neg(mat->row[var->index][off + j]), return -1); 1615 recheck = isl_tab_kill_col(tab, j); 1616 if (recheck < 0) 1617 return -1; 1618 if (recheck) 1619 --j; 1620 } 1621 if (isl_tab_mark_redundant(tab, var->index) < 0) 1622 return -1; 1623 if (tab_is_manifestly_empty(tab) && isl_tab_mark_empty(tab) < 0) 1624 return -1; 1625 return 0; 1626} 1627 1628/* Add a constraint to the tableau and allocate a row for it. 1629 * Return the index into the constraint array "con". 1630 */ 1631int isl_tab_allocate_con(struct isl_tab *tab) 1632{ 1633 int r; 1634 1635 isl_assert(tab->mat->ctx, tab->n_row < tab->mat->n_row, return -1); 1636 isl_assert(tab->mat->ctx, tab->n_con < tab->max_con, return -1); 1637 1638 r = tab->n_con; 1639 tab->con[r].index = tab->n_row; 1640 tab->con[r].is_row = 1; 1641 tab->con[r].is_nonneg = 0; 1642 tab->con[r].is_zero = 0; 1643 tab->con[r].is_redundant = 0; 1644 tab->con[r].frozen = 0; 1645 tab->con[r].negated = 0; 1646 tab->row_var[tab->n_row] = ~r; 1647 1648 tab->n_row++; 1649 tab->n_con++; 1650 if (isl_tab_push_var(tab, isl_tab_undo_allocate, &tab->con[r]) < 0) 1651 return -1; 1652 1653 return r; 1654} 1655 1656/* Add a variable to the tableau and allocate a column for it. 1657 * Return the index into the variable array "var". 1658 */ 1659int isl_tab_allocate_var(struct isl_tab *tab) 1660{ 1661 int r; 1662 int i; 1663 unsigned off = 2 + tab->M; 1664 1665 isl_assert(tab->mat->ctx, tab->n_col < tab->mat->n_col, return -1); 1666 isl_assert(tab->mat->ctx, tab->n_var < tab->max_var, return -1); 1667 1668 r = tab->n_var; 1669 tab->var[r].index = tab->n_col; 1670 tab->var[r].is_row = 0; 1671 tab->var[r].is_nonneg = 0; 1672 tab->var[r].is_zero = 0; 1673 tab->var[r].is_redundant = 0; 1674 tab->var[r].frozen = 0; 1675 tab->var[r].negated = 0; 1676 tab->col_var[tab->n_col] = r; 1677 1678 for (i = 0; i < tab->n_row; ++i) 1679 isl_int_set_si(tab->mat->row[i][off + tab->n_col], 0); 1680 1681 tab->n_var++; 1682 tab->n_col++; 1683 if (isl_tab_push_var(tab, isl_tab_undo_allocate, &tab->var[r]) < 0) 1684 return -1; 1685 1686 return r; 1687} 1688 1689/* Add a row to the tableau. The row is given as an affine combination 1690 * of the original variables and needs to be expressed in terms of the 1691 * column variables. 1692 * 1693 * We add each term in turn. 1694 * If r = n/d_r is the current sum and we need to add k x, then 1695 * if x is a column variable, we increase the numerator of 1696 * this column by k d_r 1697 * if x = f/d_x is a row variable, then the new representation of r is 1698 * 1699 * n k f d_x/g n + d_r/g k f m/d_r n + m/d_g k f 1700 * --- + --- = ------------------- = ------------------- 1701 * d_r d_r d_r d_x/g m 1702 * 1703 * with g the gcd of d_r and d_x and m the lcm of d_r and d_x. 1704 * 1705 * If tab->M is set, then, internally, each variable x is represented 1706 * as x' - M. We then also need no subtract k d_r from the coefficient of M. 1707 */ 1708int isl_tab_add_row(struct isl_tab *tab, isl_int *line) 1709{ 1710 int i; 1711 int r; 1712 isl_int *row; 1713 isl_int a, b; 1714 unsigned off = 2 + tab->M; 1715 1716 r = isl_tab_allocate_con(tab); 1717 if (r < 0) 1718 return -1; 1719 1720 isl_int_init(a); 1721 isl_int_init(b); 1722 row = tab->mat->row[tab->con[r].index]; 1723 isl_int_set_si(row[0], 1); 1724 isl_int_set(row[1], line[0]); 1725 isl_seq_clr(row + 2, tab->M + tab->n_col); 1726 for (i = 0; i < tab->n_var; ++i) { 1727 if (tab->var[i].is_zero) 1728 continue; 1729 if (tab->var[i].is_row) { 1730 isl_int_lcm(a, 1731 row[0], tab->mat->row[tab->var[i].index][0]); 1732 isl_int_swap(a, row[0]); 1733 isl_int_divexact(a, row[0], a); 1734 isl_int_divexact(b, 1735 row[0], tab->mat->row[tab->var[i].index][0]); 1736 isl_int_mul(b, b, line[1 + i]); 1737 isl_seq_combine(row + 1, a, row + 1, 1738 b, tab->mat->row[tab->var[i].index] + 1, 1739 1 + tab->M + tab->n_col); 1740 } else 1741 isl_int_addmul(row[off + tab->var[i].index], 1742 line[1 + i], row[0]); 1743 if (tab->M && i >= tab->n_param && i < tab->n_var - tab->n_div) 1744 isl_int_submul(row[2], line[1 + i], row[0]); 1745 } 1746 isl_seq_normalize(tab->mat->ctx, row, off + tab->n_col); 1747 isl_int_clear(a); 1748 isl_int_clear(b); 1749 1750 if (tab->row_sign) 1751 tab->row_sign[tab->con[r].index] = isl_tab_row_unknown; 1752 1753 return r; 1754} 1755 1756static int drop_row(struct isl_tab *tab, int row) 1757{ 1758 isl_assert(tab->mat->ctx, ~tab->row_var[row] == tab->n_con - 1, return -1); 1759 if (row != tab->n_row - 1) 1760 swap_rows(tab, row, tab->n_row - 1); 1761 tab->n_row--; 1762 tab->n_con--; 1763 return 0; 1764} 1765 1766static int drop_col(struct isl_tab *tab, int col) 1767{ 1768 isl_assert(tab->mat->ctx, tab->col_var[col] == tab->n_var - 1, return -1); 1769 if (col != tab->n_col - 1) 1770 swap_cols(tab, col, tab->n_col - 1); 1771 tab->n_col--; 1772 tab->n_var--; 1773 return 0; 1774} 1775 1776/* Add inequality "ineq" and check if it conflicts with the 1777 * previously added constraints or if it is obviously redundant. 1778 */ 1779int isl_tab_add_ineq(struct isl_tab *tab, isl_int *ineq) 1780{ 1781 int r; 1782 int sgn; 1783 isl_int cst; 1784 1785 if (!tab) 1786 return -1; 1787 if (tab->bmap) { 1788 struct isl_basic_map *bmap = tab->bmap; 1789 1790 isl_assert(tab->mat->ctx, tab->n_eq == bmap->n_eq, return -1); 1791 isl_assert(tab->mat->ctx, 1792 tab->n_con == bmap->n_eq + bmap->n_ineq, return -1); 1793 tab->bmap = isl_basic_map_add_ineq(tab->bmap, ineq); 1794 if (isl_tab_push(tab, isl_tab_undo_bmap_ineq) < 0) 1795 return -1; 1796 if (!tab->bmap) 1797 return -1; 1798 } 1799 if (tab->cone) { 1800 isl_int_init(cst); 1801 isl_int_swap(ineq[0], cst); 1802 } 1803 r = isl_tab_add_row(tab, ineq); 1804 if (tab->cone) { 1805 isl_int_swap(ineq[0], cst); 1806 isl_int_clear(cst); 1807 } 1808 if (r < 0) 1809 return -1; 1810 tab->con[r].is_nonneg = 1; 1811 if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r]) < 0) 1812 return -1; 1813 if (isl_tab_row_is_redundant(tab, tab->con[r].index)) { 1814 if (isl_tab_mark_redundant(tab, tab->con[r].index) < 0) 1815 return -1; 1816 return 0; 1817 } 1818 1819 sgn = restore_row(tab, &tab->con[r]); 1820 if (sgn < -1) 1821 return -1; 1822 if (sgn < 0) 1823 return isl_tab_mark_empty(tab); 1824 if (tab->con[r].is_row && isl_tab_row_is_redundant(tab, tab->con[r].index)) 1825 if (isl_tab_mark_redundant(tab, tab->con[r].index) < 0) 1826 return -1; 1827 return 0; 1828} 1829 1830/* Pivot a non-negative variable down until it reaches the value zero 1831 * and then pivot the variable into a column position. 1832 */ 1833static int to_col(struct isl_tab *tab, struct isl_tab_var *var) WARN_UNUSED; 1834static int to_col(struct isl_tab *tab, struct isl_tab_var *var) 1835{ 1836 int i; 1837 int row, col; 1838 unsigned off = 2 + tab->M; 1839 1840 if (!var->is_row) 1841 return 0; 1842 1843 while (isl_int_is_pos(tab->mat->row[var->index][1])) { 1844 find_pivot(tab, var, NULL, -1, &row, &col); 1845 isl_assert(tab->mat->ctx, row != -1, return -1); 1846 if (isl_tab_pivot(tab, row, col) < 0) 1847 return -1; 1848 if (!var->is_row) 1849 return 0; 1850 } 1851 1852 for (i = tab->n_dead; i < tab->n_col; ++i) 1853 if (!isl_int_is_zero(tab->mat->row[var->index][off + i])) 1854 break; 1855 1856 isl_assert(tab->mat->ctx, i < tab->n_col, return -1); 1857 if (isl_tab_pivot(tab, var->index, i) < 0) 1858 return -1; 1859 1860 return 0; 1861} 1862 1863/* We assume Gaussian elimination has been performed on the equalities. 1864 * The equalities can therefore never conflict. 1865 * Adding the equalities is currently only really useful for a later call 1866 * to isl_tab_ineq_type. 1867 */ 1868static struct isl_tab *add_eq(struct isl_tab *tab, isl_int *eq) 1869{ 1870 int i; 1871 int r; 1872 1873 if (!tab) 1874 return NULL; 1875 r = isl_tab_add_row(tab, eq); 1876 if (r < 0) 1877 goto error; 1878 1879 r = tab->con[r].index; 1880 i = isl_seq_first_non_zero(tab->mat->row[r] + 2 + tab->M + tab->n_dead, 1881 tab->n_col - tab->n_dead); 1882 isl_assert(tab->mat->ctx, i >= 0, goto error); 1883 i += tab->n_dead; 1884 if (isl_tab_pivot(tab, r, i) < 0) 1885 goto error; 1886 if (isl_tab_kill_col(tab, i) < 0) 1887 goto error; 1888 tab->n_eq++; 1889 1890 return tab; 1891error: 1892 isl_tab_free(tab); 1893 return NULL; 1894} 1895 1896static int row_is_manifestly_zero(struct isl_tab *tab, int row) 1897{ 1898 unsigned off = 2 + tab->M; 1899 1900 if (!isl_int_is_zero(tab->mat->row[row][1])) 1901 return 0; 1902 if (tab->M && !isl_int_is_zero(tab->mat->row[row][2])) 1903 return 0; 1904 return isl_seq_first_non_zero(tab->mat->row[row] + off + tab->n_dead, 1905 tab->n_col - tab->n_dead) == -1; 1906} 1907 1908/* Add an equality that is known to be valid for the given tableau. 1909 */ 1910int isl_tab_add_valid_eq(struct isl_tab *tab, isl_int *eq) 1911{ 1912 struct isl_tab_var *var; 1913 int r; 1914 1915 if (!tab) 1916 return -1; 1917 r = isl_tab_add_row(tab, eq); 1918 if (r < 0) 1919 return -1; 1920 1921 var = &tab->con[r]; 1922 r = var->index; 1923 if (row_is_manifestly_zero(tab, r)) { 1924 var->is_zero = 1; 1925 if (isl_tab_mark_redundant(tab, r) < 0) 1926 return -1; 1927 return 0; 1928 } 1929 1930 if (isl_int_is_neg(tab->mat->row[r][1])) { 1931 isl_seq_neg(tab->mat->row[r] + 1, tab->mat->row[r] + 1, 1932 1 + tab->n_col); 1933 var->negated = 1; 1934 } 1935 var->is_nonneg = 1; 1936 if (to_col(tab, var) < 0) 1937 return -1; 1938 var->is_nonneg = 0; 1939 if (isl_tab_kill_col(tab, var->index) < 0) 1940 return -1; 1941 1942 return 0; 1943} 1944 1945static int add_zero_row(struct isl_tab *tab) 1946{ 1947 int r; 1948 isl_int *row; 1949 1950 r = isl_tab_allocate_con(tab); 1951 if (r < 0) 1952 return -1; 1953 1954 row = tab->mat->row[tab->con[r].index]; 1955 isl_seq_clr(row + 1, 1 + tab->M + tab->n_col); 1956 isl_int_set_si(row[0], 1); 1957 1958 return r; 1959} 1960 1961/* Add equality "eq" and check if it conflicts with the 1962 * previously added constraints or if it is obviously redundant. 1963 */ 1964int isl_tab_add_eq(struct isl_tab *tab, isl_int *eq) 1965{ 1966 struct isl_tab_undo *snap = NULL; 1967 struct isl_tab_var *var; 1968 int r; 1969 int row; 1970 int sgn; 1971 isl_int cst; 1972 1973 if (!tab) 1974 return -1; 1975 isl_assert(tab->mat->ctx, !tab->M, return -1); 1976 1977 if (tab->need_undo) 1978 snap = isl_tab_snap(tab); 1979 1980 if (tab->cone) { 1981 isl_int_init(cst); 1982 isl_int_swap(eq[0], cst); 1983 } 1984 r = isl_tab_add_row(tab, eq); 1985 if (tab->cone) { 1986 isl_int_swap(eq[0], cst); 1987 isl_int_clear(cst); 1988 } 1989 if (r < 0) 1990 return -1; 1991 1992 var = &tab->con[r]; 1993 row = var->index; 1994 if (row_is_manifestly_zero(tab, row)) { 1995 if (snap) { 1996 if (isl_tab_rollback(tab, snap) < 0) 1997 return -1; 1998 } else 1999 drop_row(tab, row); 2000 return 0; 2001 } 2002 2003 if (tab->bmap) { 2004 tab->bmap = isl_basic_map_add_ineq(tab->bmap, eq); 2005 if (isl_tab_push(tab, isl_tab_undo_bmap_ineq) < 0) 2006 return -1; 2007 isl_seq_neg(eq, eq, 1 + tab->n_var); 2008 tab->bmap = isl_basic_map_add_ineq(tab->bmap, eq); 2009 isl_seq_neg(eq, eq, 1 + tab->n_var); 2010 if (isl_tab_push(tab, isl_tab_undo_bmap_ineq) < 0) 2011 return -1; 2012 if (!tab->bmap) 2013 return -1; 2014 if (add_zero_row(tab) < 0) 2015 return -1; 2016 } 2017 2018 sgn = isl_int_sgn(tab->mat->row[row][1]); 2019 2020 if (sgn > 0) { 2021 isl_seq_neg(tab->mat->row[row] + 1, tab->mat->row[row] + 1, 2022 1 + tab->n_col); 2023 var->negated = 1; 2024 sgn = -1; 2025 } 2026 2027 if (sgn < 0) { 2028 sgn = sign_of_max(tab, var); 2029 if (sgn < -1) 2030 return -1; 2031 if (sgn < 0) { 2032 if (isl_tab_mark_empty(tab) < 0) 2033 return -1; 2034 return 0; 2035 } 2036 } 2037 2038 var->is_nonneg = 1; 2039 if (to_col(tab, var) < 0) 2040 return -1; 2041 var->is_nonneg = 0; 2042 if (isl_tab_kill_col(tab, var->index) < 0) 2043 return -1; 2044 2045 return 0; 2046} 2047 2048/* Construct and return an inequality that expresses an upper bound 2049 * on the given div. 2050 * In particular, if the div is given by 2051 * 2052 * d = floor(e/m) 2053 * 2054 * then the inequality expresses 2055 * 2056 * m d <= e 2057 */ 2058static struct isl_vec *ineq_for_div(struct isl_basic_map *bmap, unsigned div) 2059{ 2060 unsigned total; 2061 unsigned div_pos; 2062 struct isl_vec *ineq; 2063 2064 if (!bmap) 2065 return NULL; 2066 2067 total = isl_basic_map_total_dim(bmap); 2068 div_pos = 1 + total - bmap->n_div + div; 2069 2070 ineq = isl_vec_alloc(bmap->ctx, 1 + total); 2071 if (!ineq) 2072 return NULL; 2073 2074 isl_seq_cpy(ineq->el, bmap->div[div] + 1, 1 + total); 2075 isl_int_neg(ineq->el[div_pos], bmap->div[div][0]); 2076 return ineq; 2077} 2078 2079/* For a div d = floor(f/m), add the constraints 2080 * 2081 * f - m d >= 0 2082 * -(f-(m-1)) + m d >= 0 2083 * 2084 * Note that the second constraint is the negation of 2085 * 2086 * f - m d >= m 2087 * 2088 * If add_ineq is not NULL, then this function is used 2089 * instead of isl_tab_add_ineq to effectively add the inequalities. 2090 */ 2091static int add_div_constraints(struct isl_tab *tab, unsigned div, 2092 int (*add_ineq)(void *user, isl_int *), void *user) 2093{ 2094 unsigned total; 2095 unsigned div_pos; 2096 struct isl_vec *ineq; 2097 2098 total = isl_basic_map_total_dim(tab->bmap); 2099 div_pos = 1 + total - tab->bmap->n_div + div; 2100 2101 ineq = ineq_for_div(tab->bmap, div); 2102 if (!ineq) 2103 goto error; 2104 2105 if (add_ineq) { 2106 if (add_ineq(user, ineq->el) < 0) 2107 goto error; 2108 } else { 2109 if (isl_tab_add_ineq(tab, ineq->el) < 0) 2110 goto error; 2111 } 2112 2113 isl_seq_neg(ineq->el, tab->bmap->div[div] + 1, 1 + total); 2114 isl_int_set(ineq->el[div_pos], tab->bmap->div[div][0]); 2115 isl_int_add(ineq->el[0], ineq->el[0], ineq->el[div_pos]); 2116 isl_int_sub_ui(ineq->el[0], ineq->el[0], 1); 2117 2118 if (add_ineq) { 2119 if (add_ineq(user, ineq->el) < 0) 2120 goto error; 2121 } else { 2122 if (isl_tab_add_ineq(tab, ineq->el) < 0) 2123 goto error; 2124 } 2125 2126 isl_vec_free(ineq); 2127 2128 return 0; 2129error: 2130 isl_vec_free(ineq); 2131 return -1; 2132} 2133 2134/* Check whether the div described by "div" is obviously non-negative. 2135 * If we are using a big parameter, then we will encode the div 2136 * as div' = M + div, which is always non-negative. 2137 * Otherwise, we check whether div is a non-negative affine combination 2138 * of non-negative variables. 2139 */ 2140static int div_is_nonneg(struct isl_tab *tab, __isl_keep isl_vec *div) 2141{ 2142 int i; 2143 2144 if (tab->M) 2145 return 1; 2146 2147 if (isl_int_is_neg(div->el[1])) 2148 return 0; 2149 2150 for (i = 0; i < tab->n_var; ++i) { 2151 if (isl_int_is_neg(div->el[2 + i])) 2152 return 0; 2153 if (isl_int_is_zero(div->el[2 + i])) 2154 continue; 2155 if (!tab->var[i].is_nonneg) 2156 return 0; 2157 } 2158 2159 return 1; 2160} 2161 2162/* Add an extra div, prescribed by "div" to the tableau and 2163 * the associated bmap (which is assumed to be non-NULL). 2164 * 2165 * If add_ineq is not NULL, then this function is used instead 2166 * of isl_tab_add_ineq to add the div constraints. 2167 * This complication is needed because the code in isl_tab_pip 2168 * wants to perform some extra processing when an inequality 2169 * is added to the tableau. 2170 */ 2171int isl_tab_add_div(struct isl_tab *tab, __isl_keep isl_vec *div, 2172 int (*add_ineq)(void *user, isl_int *), void *user) 2173{ 2174 int r; 2175 int k; 2176 int nonneg; 2177 2178 if (!tab || !div) 2179 return -1; 2180 2181 isl_assert(tab->mat->ctx, tab->bmap, return -1); 2182 2183 nonneg = div_is_nonneg(tab, div); 2184 2185 if (isl_tab_extend_cons(tab, 3) < 0) 2186 return -1; 2187 if (isl_tab_extend_vars(tab, 1) < 0) 2188 return -1; 2189 r = isl_tab_allocate_var(tab); 2190 if (r < 0) 2191 return -1; 2192 2193 if (nonneg) 2194 tab->var[r].is_nonneg = 1; 2195 2196 tab->bmap = isl_basic_map_extend_space(tab->bmap, 2197 isl_basic_map_get_space(tab->bmap), 1, 0, 2); 2198 k = isl_basic_map_alloc_div(tab->bmap); 2199 if (k < 0) 2200 return -1; 2201 isl_seq_cpy(tab->bmap->div[k], div->el, div->size); 2202 if (isl_tab_push(tab, isl_tab_undo_bmap_div) < 0) 2203 return -1; 2204 2205 if (add_div_constraints(tab, k, add_ineq, user) < 0) 2206 return -1; 2207 2208 return r; 2209} 2210 2211/* If "track" is set, then we want to keep track of all constraints in tab 2212 * in its bmap field. This field is initialized from a copy of "bmap", 2213 * so we need to make sure that all constraints in "bmap" also appear 2214 * in the constructed tab. 2215 */ 2216__isl_give struct isl_tab *isl_tab_from_basic_map( 2217 __isl_keep isl_basic_map *bmap, int track) 2218{ 2219 int i; 2220 struct isl_tab *tab; 2221 2222 if (!bmap) 2223 return NULL; 2224 tab = isl_tab_alloc(bmap->ctx, 2225 isl_basic_map_total_dim(bmap) + bmap->n_ineq + 1, 2226 isl_basic_map_total_dim(bmap), 0); 2227 if (!tab) 2228 return NULL; 2229 tab->preserve = track; 2230 tab->rational = ISL_F_ISSET(bmap, ISL_BASIC_MAP_RATIONAL); 2231 if (ISL_F_ISSET(bmap, ISL_BASIC_MAP_EMPTY)) { 2232 if (isl_tab_mark_empty(tab) < 0) 2233 goto error; 2234 goto done; 2235 } 2236 for (i = 0; i < bmap->n_eq; ++i) { 2237 tab = add_eq(tab, bmap->eq[i]); 2238 if (!tab) 2239 return tab; 2240 } 2241 for (i = 0; i < bmap->n_ineq; ++i) { 2242 if (isl_tab_add_ineq(tab, bmap->ineq[i]) < 0) 2243 goto error; 2244 if (tab->empty) 2245 goto done; 2246 } 2247done: 2248 if (track && isl_tab_track_bmap(tab, isl_basic_map_copy(bmap)) < 0) 2249 goto error; 2250 return tab; 2251error: 2252 isl_tab_free(tab); 2253 return NULL; 2254} 2255 2256__isl_give struct isl_tab *isl_tab_from_basic_set( 2257 __isl_keep isl_basic_set *bset, int track) 2258{ 2259 return isl_tab_from_basic_map(bset, track); 2260} 2261 2262/* Construct a tableau corresponding to the recession cone of "bset". 2263 */ 2264struct isl_tab *isl_tab_from_recession_cone(__isl_keep isl_basic_set *bset, 2265 int parametric) 2266{ 2267 isl_int cst; 2268 int i; 2269 struct isl_tab *tab; 2270 unsigned offset = 0; 2271 2272 if (!bset) 2273 return NULL; 2274 if (parametric) 2275 offset = isl_basic_set_dim(bset, isl_dim_param); 2276 tab = isl_tab_alloc(bset->ctx, bset->n_eq + bset->n_ineq, 2277 isl_basic_set_total_dim(bset) - offset, 0); 2278 if (!tab) 2279 return NULL; 2280 tab->rational = ISL_F_ISSET(bset, ISL_BASIC_SET_RATIONAL); 2281 tab->cone = 1; 2282 2283 isl_int_init(cst); 2284 for (i = 0; i < bset->n_eq; ++i) { 2285 isl_int_swap(bset->eq[i][offset], cst); 2286 if (offset > 0) { 2287 if (isl_tab_add_eq(tab, bset->eq[i] + offset) < 0) 2288 goto error; 2289 } else 2290 tab = add_eq(tab, bset->eq[i]); 2291 isl_int_swap(bset->eq[i][offset], cst); 2292 if (!tab) 2293 goto done; 2294 } 2295 for (i = 0; i < bset->n_ineq; ++i) { 2296 int r; 2297 isl_int_swap(bset->ineq[i][offset], cst); 2298 r = isl_tab_add_row(tab, bset->ineq[i] + offset); 2299 isl_int_swap(bset->ineq[i][offset], cst); 2300 if (r < 0) 2301 goto error; 2302 tab->con[r].is_nonneg = 1; 2303 if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r]) < 0) 2304 goto error; 2305 } 2306done: 2307 isl_int_clear(cst); 2308 return tab; 2309error: 2310 isl_int_clear(cst); 2311 isl_tab_free(tab); 2312 return NULL; 2313} 2314 2315/* Assuming "tab" is the tableau of a cone, check if the cone is 2316 * bounded, i.e., if it is empty or only contains the origin. 2317 */ 2318int isl_tab_cone_is_bounded(struct isl_tab *tab) 2319{ 2320 int i; 2321 2322 if (!tab) 2323 return -1; 2324 if (tab->empty) 2325 return 1; 2326 if (tab->n_dead == tab->n_col) 2327 return 1; 2328 2329 for (;;) { 2330 for (i = tab->n_redundant; i < tab->n_row; ++i) { 2331 struct isl_tab_var *var; 2332 int sgn; 2333 var = isl_tab_var_from_row(tab, i); 2334 if (!var->is_nonneg) 2335 continue; 2336 sgn = sign_of_max(tab, var); 2337 if (sgn < -1) 2338 return -1; 2339 if (sgn != 0) 2340 return 0; 2341 if (close_row(tab, var) < 0) 2342 return -1; 2343 break; 2344 } 2345 if (tab->n_dead == tab->n_col) 2346 return 1; 2347 if (i == tab->n_row) 2348 return 0; 2349 } 2350} 2351 2352int isl_tab_sample_is_integer(struct isl_tab *tab) 2353{ 2354 int i; 2355 2356 if (!tab) 2357 return -1; 2358 2359 for (i = 0; i < tab->n_var; ++i) { 2360 int row; 2361 if (!tab->var[i].is_row) 2362 continue; 2363 row = tab->var[i].index; 2364 if (!isl_int_is_divisible_by(tab->mat->row[row][1], 2365 tab->mat->row[row][0])) 2366 return 0; 2367 } 2368 return 1; 2369} 2370 2371static struct isl_vec *extract_integer_sample(struct isl_tab *tab) 2372{ 2373 int i; 2374 struct isl_vec *vec; 2375 2376 vec = isl_vec_alloc(tab->mat->ctx, 1 + tab->n_var); 2377 if (!vec) 2378 return NULL; 2379 2380 isl_int_set_si(vec->block.data[0], 1); 2381 for (i = 0; i < tab->n_var; ++i) { 2382 if (!tab->var[i].is_row) 2383 isl_int_set_si(vec->block.data[1 + i], 0); 2384 else { 2385 int row = tab->var[i].index; 2386 isl_int_divexact(vec->block.data[1 + i], 2387 tab->mat->row[row][1], tab->mat->row[row][0]); 2388 } 2389 } 2390 2391 return vec; 2392} 2393 2394struct isl_vec *isl_tab_get_sample_value(struct isl_tab *tab) 2395{ 2396 int i; 2397 struct isl_vec *vec; 2398 isl_int m; 2399 2400 if (!tab) 2401 return NULL; 2402 2403 vec = isl_vec_alloc(tab->mat->ctx, 1 + tab->n_var); 2404 if (!vec) 2405 return NULL; 2406 2407 isl_int_init(m); 2408 2409 isl_int_set_si(vec->block.data[0], 1); 2410 for (i = 0; i < tab->n_var; ++i) { 2411 int row; 2412 if (!tab->var[i].is_row) { 2413 isl_int_set_si(vec->block.data[1 + i], 0); 2414 continue; 2415 } 2416 row = tab->var[i].index; 2417 isl_int_gcd(m, vec->block.data[0], tab->mat->row[row][0]); 2418 isl_int_divexact(m, tab->mat->row[row][0], m); 2419 isl_seq_scale(vec->block.data, vec->block.data, m, 1 + i); 2420 isl_int_divexact(m, vec->block.data[0], tab->mat->row[row][0]); 2421 isl_int_mul(vec->block.data[1 + i], m, tab->mat->row[row][1]); 2422 } 2423 vec = isl_vec_normalize(vec); 2424 2425 isl_int_clear(m); 2426 return vec; 2427} 2428 2429/* Update "bmap" based on the results of the tableau "tab". 2430 * In particular, implicit equalities are made explicit, redundant constraints 2431 * are removed and if the sample value happens to be integer, it is stored 2432 * in "bmap" (unless "bmap" already had an integer sample). 2433 * 2434 * The tableau is assumed to have been created from "bmap" using 2435 * isl_tab_from_basic_map. 2436 */ 2437struct isl_basic_map *isl_basic_map_update_from_tab(struct isl_basic_map *bmap, 2438 struct isl_tab *tab) 2439{ 2440 int i; 2441 unsigned n_eq; 2442 2443 if (!bmap) 2444 return NULL; 2445 if (!tab) 2446 return bmap; 2447 2448 n_eq = tab->n_eq; 2449 if (tab->empty) 2450 bmap = isl_basic_map_set_to_empty(bmap); 2451 else 2452 for (i = bmap->n_ineq - 1; i >= 0; --i) { 2453 if (isl_tab_is_equality(tab, n_eq + i)) 2454 isl_basic_map_inequality_to_equality(bmap, i); 2455 else if (isl_tab_is_redundant(tab, n_eq + i)) 2456 isl_basic_map_drop_inequality(bmap, i); 2457 } 2458 if (bmap->n_eq != n_eq) 2459 isl_basic_map_gauss(bmap, NULL); 2460 if (!tab->rational && 2461 !bmap->sample && isl_tab_sample_is_integer(tab)) 2462 bmap->sample = extract_integer_sample(tab); 2463 return bmap; 2464} 2465 2466struct isl_basic_set *isl_basic_set_update_from_tab(struct isl_basic_set *bset, 2467 struct isl_tab *tab) 2468{ 2469 return (struct isl_basic_set *)isl_basic_map_update_from_tab( 2470 (struct isl_basic_map *)bset, tab); 2471} 2472 2473/* Given a non-negative variable "var", add a new non-negative variable 2474 * that is the opposite of "var", ensuring that var can only attain the 2475 * value zero. 2476 * If var = n/d is a row variable, then the new variable = -n/d. 2477 * If var is a column variables, then the new variable = -var. 2478 * If the new variable cannot attain non-negative values, then 2479 * the resulting tableau is empty. 2480 * Otherwise, we know the value will be zero and we close the row. 2481 */ 2482static int cut_to_hyperplane(struct isl_tab *tab, struct isl_tab_var *var) 2483{ 2484 unsigned r; 2485 isl_int *row; 2486 int sgn; 2487 unsigned off = 2 + tab->M; 2488 2489 if (var->is_zero) 2490 return 0; 2491 isl_assert(tab->mat->ctx, !var->is_redundant, return -1); 2492 isl_assert(tab->mat->ctx, var->is_nonneg, return -1); 2493 2494 if (isl_tab_extend_cons(tab, 1) < 0) 2495 return -1; 2496 2497 r = tab->n_con; 2498 tab->con[r].index = tab->n_row; 2499 tab->con[r].is_row = 1; 2500 tab->con[r].is_nonneg = 0; 2501 tab->con[r].is_zero = 0; 2502 tab->con[r].is_redundant = 0; 2503 tab->con[r].frozen = 0; 2504 tab->con[r].negated = 0; 2505 tab->row_var[tab->n_row] = ~r; 2506 row = tab->mat->row[tab->n_row]; 2507 2508 if (var->is_row) { 2509 isl_int_set(row[0], tab->mat->row[var->index][0]); 2510 isl_seq_neg(row + 1, 2511 tab->mat->row[var->index] + 1, 1 + tab->n_col); 2512 } else { 2513 isl_int_set_si(row[0], 1); 2514 isl_seq_clr(row + 1, 1 + tab->n_col); 2515 isl_int_set_si(row[off + var->index], -1); 2516 } 2517 2518 tab->n_row++; 2519 tab->n_con++; 2520 if (isl_tab_push_var(tab, isl_tab_undo_allocate, &tab->con[r]) < 0) 2521 return -1; 2522 2523 sgn = sign_of_max(tab, &tab->con[r]); 2524 if (sgn < -1) 2525 return -1; 2526 if (sgn < 0) { 2527 if (isl_tab_mark_empty(tab) < 0) 2528 return -1; 2529 return 0; 2530 } 2531 tab->con[r].is_nonneg = 1; 2532 if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r]) < 0) 2533 return -1; 2534 /* sgn == 0 */ 2535 if (close_row(tab, &tab->con[r]) < 0) 2536 return -1; 2537 2538 return 0; 2539} 2540 2541/* Given a tableau "tab" and an inequality constraint "con" of the tableau, 2542 * relax the inequality by one. That is, the inequality r >= 0 is replaced 2543 * by r' = r + 1 >= 0. 2544 * If r is a row variable, we simply increase the constant term by one 2545 * (taking into account the denominator). 2546 * If r is a column variable, then we need to modify each row that 2547 * refers to r = r' - 1 by substituting this equality, effectively 2548 * subtracting the coefficient of the column from the constant. 2549 * We should only do this if the minimum is manifestly unbounded, 2550 * however. Otherwise, we may end up with negative sample values 2551 * for non-negative variables. 2552 * So, if r is a column variable with a minimum that is not 2553 * manifestly unbounded, then we need to move it to a row. 2554 * However, the sample value of this row may be negative, 2555 * even after the relaxation, so we need to restore it. 2556 * We therefore prefer to pivot a column up to a row, if possible. 2557 */ 2558struct isl_tab *isl_tab_relax(struct isl_tab *tab, int con) 2559{ 2560 struct isl_tab_var *var; 2561 unsigned off = 2 + tab->M; 2562 2563 if (!tab) 2564 return NULL; 2565 2566 var = &tab->con[con]; 2567 2568 if (var->is_row && (var->index < 0 || var->index < tab->n_redundant)) 2569 isl_die(tab->mat->ctx, isl_error_invalid, 2570 "cannot relax redundant constraint", goto error); 2571 if (!var->is_row && (var->index < 0 || var->index < tab->n_dead)) 2572 isl_die(tab->mat->ctx, isl_error_invalid, 2573 "cannot relax dead constraint", goto error); 2574 2575 if (!var->is_row && !max_is_manifestly_unbounded(tab, var)) 2576 if (to_row(tab, var, 1) < 0) 2577 goto error; 2578 if (!var->is_row && !min_is_manifestly_unbounded(tab, var)) 2579 if (to_row(tab, var, -1) < 0) 2580 goto error; 2581 2582 if (var->is_row) { 2583 isl_int_add(tab->mat->row[var->index][1], 2584 tab->mat->row[var->index][1], tab->mat->row[var->index][0]); 2585 if (restore_row(tab, var) < 0) 2586 goto error; 2587 } else { 2588 int i; 2589 2590 for (i = 0; i < tab->n_row; ++i) { 2591 if (isl_int_is_zero(tab->mat->row[i][off + var->index])) 2592 continue; 2593 isl_int_sub(tab->mat->row[i][1], tab->mat->row[i][1], 2594 tab->mat->row[i][off + var->index]); 2595 } 2596 2597 } 2598 2599 if (isl_tab_push_var(tab, isl_tab_undo_relax, var) < 0) 2600 goto error; 2601 2602 return tab; 2603error: 2604 isl_tab_free(tab); 2605 return NULL; 2606} 2607 2608/* Remove the sign constraint from constraint "con". 2609 * 2610 * If the constraint variable was originally marked non-negative, 2611 * then we make sure we mark it non-negative again during rollback. 2612 */ 2613int isl_tab_unrestrict(struct isl_tab *tab, int con) 2614{ 2615 struct isl_tab_var *var; 2616 2617 if (!tab) 2618 return -1; 2619 2620 var = &tab->con[con]; 2621 if (!var->is_nonneg) 2622 return 0; 2623 2624 var->is_nonneg = 0; 2625 if (isl_tab_push_var(tab, isl_tab_undo_unrestrict, var) < 0) 2626 return -1; 2627 2628 return 0; 2629} 2630 2631int isl_tab_select_facet(struct isl_tab *tab, int con) 2632{ 2633 if (!tab) 2634 return -1; 2635 2636 return cut_to_hyperplane(tab, &tab->con[con]); 2637} 2638 2639static int may_be_equality(struct isl_tab *tab, int row) 2640{ 2641 return tab->rational ? isl_int_is_zero(tab->mat->row[row][1]) 2642 : isl_int_lt(tab->mat->row[row][1], 2643 tab->mat->row[row][0]); 2644} 2645 2646/* Check for (near) equalities among the constraints. 2647 * A constraint is an equality if it is non-negative and if 2648 * its maximal value is either 2649 * - zero (in case of rational tableaus), or 2650 * - strictly less than 1 (in case of integer tableaus) 2651 * 2652 * We first mark all non-redundant and non-dead variables that 2653 * are not frozen and not obviously not an equality. 2654 * Then we iterate over all marked variables if they can attain 2655 * any values larger than zero or at least one. 2656 * If the maximal value is zero, we mark any column variables 2657 * that appear in the row as being zero and mark the row as being redundant. 2658 * Otherwise, if the maximal value is strictly less than one (and the 2659 * tableau is integer), then we restrict the value to being zero 2660 * by adding an opposite non-negative variable. 2661 */ 2662int isl_tab_detect_implicit_equalities(struct isl_tab *tab) 2663{ 2664 int i; 2665 unsigned n_marked; 2666 2667 if (!tab) 2668 return -1; 2669 if (tab->empty) 2670 return 0; 2671 if (tab->n_dead == tab->n_col) 2672 return 0; 2673 2674 n_marked = 0; 2675 for (i = tab->n_redundant; i < tab->n_row; ++i) { 2676 struct isl_tab_var *var = isl_tab_var_from_row(tab, i); 2677 var->marked = !var->frozen && var->is_nonneg && 2678 may_be_equality(tab, i); 2679 if (var->marked) 2680 n_marked++; 2681 } 2682 for (i = tab->n_dead; i < tab->n_col; ++i) { 2683 struct isl_tab_var *var = var_from_col(tab, i); 2684 var->marked = !var->frozen && var->is_nonneg; 2685 if (var->marked) 2686 n_marked++; 2687 } 2688 while (n_marked) { 2689 struct isl_tab_var *var; 2690 int sgn; 2691 for (i = tab->n_redundant; i < tab->n_row; ++i) { 2692 var = isl_tab_var_from_row(tab, i); 2693 if (var->marked) 2694 break; 2695 } 2696 if (i == tab->n_row) { 2697 for (i = tab->n_dead; i < tab->n_col; ++i) { 2698 var = var_from_col(tab, i); 2699 if (var->marked) 2700 break; 2701 } 2702 if (i == tab->n_col) 2703 break; 2704 } 2705 var->marked = 0; 2706 n_marked--; 2707 sgn = sign_of_max(tab, var); 2708 if (sgn < 0) 2709 return -1; 2710 if (sgn == 0) { 2711 if (close_row(tab, var) < 0) 2712 return -1; 2713 } else if (!tab->rational && !at_least_one(tab, var)) { 2714 if (cut_to_hyperplane(tab, var) < 0) 2715 return -1; 2716 return isl_tab_detect_implicit_equalities(tab); 2717 } 2718 for (i = tab->n_redundant; i < tab->n_row; ++i) { 2719 var = isl_tab_var_from_row(tab, i); 2720 if (!var->marked) 2721 continue; 2722 if (may_be_equality(tab, i)) 2723 continue; 2724 var->marked = 0; 2725 n_marked--; 2726 } 2727 } 2728 2729 return 0; 2730} 2731 2732/* Update the element of row_var or col_var that corresponds to 2733 * constraint tab->con[i] to a move from position "old" to position "i". 2734 */ 2735static int update_con_after_move(struct isl_tab *tab, int i, int old) 2736{ 2737 int *p; 2738 int index; 2739 2740 index = tab->con[i].index; 2741 if (index == -1) 2742 return 0; 2743 p = tab->con[i].is_row ? tab->row_var : tab->col_var; 2744 if (p[index] != ~old) 2745 isl_die(tab->mat->ctx, isl_error_internal, 2746 "broken internal state", return -1); 2747 p[index] = ~i; 2748 2749 return 0; 2750} 2751 2752/* Rotate the "n" constraints starting at "first" to the right, 2753 * putting the last constraint in the position of the first constraint. 2754 */ 2755static int rotate_constraints(struct isl_tab *tab, int first, int n) 2756{ 2757 int i, last; 2758 struct isl_tab_var var; 2759 2760 if (n <= 1) 2761 return 0; 2762 2763 last = first + n - 1; 2764 var = tab->con[last]; 2765 for (i = last; i > first; --i) { 2766 tab->con[i] = tab->con[i - 1]; 2767 if (update_con_after_move(tab, i, i - 1) < 0) 2768 return -1; 2769 } 2770 tab->con[first] = var; 2771 if (update_con_after_move(tab, first, last) < 0) 2772 return -1; 2773 2774 return 0; 2775} 2776 2777/* Make the equalities that are implicit in "bmap" but that have been 2778 * detected in the corresponding "tab" explicit in "bmap" and update 2779 * "tab" to reflect the new order of the constraints. 2780 * 2781 * In particular, if inequality i is an implicit equality then 2782 * isl_basic_map_inequality_to_equality will move the inequality 2783 * in front of the other equality and it will move the last inequality 2784 * in the position of inequality i. 2785 * In the tableau, the inequalities of "bmap" are stored after the equalities 2786 * and so the original order 2787 * 2788 * E E E E E A A A I B B B B L 2789 * 2790 * is changed into 2791 * 2792 * I E E E E E A A A L B B B B 2793 * 2794 * where I is the implicit equality, the E are equalities, 2795 * the A inequalities before I, the B inequalities after I and 2796 * L the last inequality. 2797 * We therefore need to rotate to the right two sets of constraints, 2798 * those up to and including I and those after I. 2799 * 2800 * If "tab" contains any constraints that are not in "bmap" then they 2801 * appear after those in "bmap" and they should be left untouched. 2802 * 2803 * Note that this function leaves "bmap" in a temporary state 2804 * as it does not call isl_basic_map_gauss. Calling this function 2805 * is the responsibility of the caller. 2806 */ 2807__isl_give isl_basic_map *isl_tab_make_equalities_explicit(struct isl_tab *tab, 2808 __isl_take isl_basic_map *bmap) 2809{ 2810 int i; 2811 2812 if (!tab || !bmap) 2813 return isl_basic_map_free(bmap); 2814 if (tab->empty) 2815 return bmap; 2816 2817 for (i = bmap->n_ineq - 1; i >= 0; --i) { 2818 if (!isl_tab_is_equality(tab, bmap->n_eq + i)) 2819 continue; 2820 isl_basic_map_inequality_to_equality(bmap, i); 2821 if (rotate_constraints(tab, 0, tab->n_eq + i + 1) < 0) 2822 return isl_basic_map_free(bmap); 2823 if (rotate_constraints(tab, tab->n_eq + i + 1, 2824 bmap->n_ineq - i) < 0) 2825 return isl_basic_map_free(bmap); 2826 tab->n_eq++; 2827 } 2828 2829 return bmap; 2830} 2831 2832static int con_is_redundant(struct isl_tab *tab, struct isl_tab_var *var) 2833{ 2834 if (!tab) 2835 return -1; 2836 if (tab->rational) { 2837 int sgn = sign_of_min(tab, var); 2838 if (sgn < -1) 2839 return -1; 2840 return sgn >= 0; 2841 } else { 2842 int irred = isl_tab_min_at_most_neg_one(tab, var); 2843 if (irred < 0) 2844 return -1; 2845 return !irred; 2846 } 2847} 2848 2849/* Check for (near) redundant constraints. 2850 * A constraint is redundant if it is non-negative and if 2851 * its minimal value (temporarily ignoring the non-negativity) is either 2852 * - zero (in case of rational tableaus), or 2853 * - strictly larger than -1 (in case of integer tableaus) 2854 * 2855 * We first mark all non-redundant and non-dead variables that 2856 * are not frozen and not obviously negatively unbounded. 2857 * Then we iterate over all marked variables if they can attain 2858 * any values smaller than zero or at most negative one. 2859 * If not, we mark the row as being redundant (assuming it hasn't 2860 * been detected as being obviously redundant in the mean time). 2861 */ 2862int isl_tab_detect_redundant(struct isl_tab *tab) 2863{ 2864 int i; 2865 unsigned n_marked; 2866 2867 if (!tab) 2868 return -1; 2869 if (tab->empty) 2870 return 0; 2871 if (tab->n_redundant == tab->n_row) 2872 return 0; 2873 2874 n_marked = 0; 2875 for (i = tab->n_redundant; i < tab->n_row; ++i) { 2876 struct isl_tab_var *var = isl_tab_var_from_row(tab, i); 2877 var->marked = !var->frozen && var->is_nonneg; 2878 if (var->marked) 2879 n_marked++; 2880 } 2881 for (i = tab->n_dead; i < tab->n_col; ++i) { 2882 struct isl_tab_var *var = var_from_col(tab, i); 2883 var->marked = !var->frozen && var->is_nonneg && 2884 !min_is_manifestly_unbounded(tab, var); 2885 if (var->marked) 2886 n_marked++; 2887 } 2888 while (n_marked) { 2889 struct isl_tab_var *var; 2890 int red; 2891 for (i = tab->n_redundant; i < tab->n_row; ++i) { 2892 var = isl_tab_var_from_row(tab, i); 2893 if (var->marked) 2894 break; 2895 } 2896 if (i == tab->n_row) { 2897 for (i = tab->n_dead; i < tab->n_col; ++i) { 2898 var = var_from_col(tab, i); 2899 if (var->marked) 2900 break; 2901 } 2902 if (i == tab->n_col) 2903 break; 2904 } 2905 var->marked = 0; 2906 n_marked--; 2907 red = con_is_redundant(tab, var); 2908 if (red < 0) 2909 return -1; 2910 if (red && !var->is_redundant) 2911 if (isl_tab_mark_redundant(tab, var->index) < 0) 2912 return -1; 2913 for (i = tab->n_dead; i < tab->n_col; ++i) { 2914 var = var_from_col(tab, i); 2915 if (!var->marked) 2916 continue; 2917 if (!min_is_manifestly_unbounded(tab, var)) 2918 continue; 2919 var->marked = 0; 2920 n_marked--; 2921 } 2922 } 2923 2924 return 0; 2925} 2926 2927int isl_tab_is_equality(struct isl_tab *tab, int con) 2928{ 2929 int row; 2930 unsigned off; 2931 2932 if (!tab) 2933 return -1; 2934 if (tab->con[con].is_zero) 2935 return 1; 2936 if (tab->con[con].is_redundant) 2937 return 0; 2938 if (!tab->con[con].is_row) 2939 return tab->con[con].index < tab->n_dead; 2940 2941 row = tab->con[con].index; 2942 2943 off = 2 + tab->M; 2944 return isl_int_is_zero(tab->mat->row[row][1]) && 2945 (!tab->M || isl_int_is_zero(tab->mat->row[row][2])) && 2946 isl_seq_first_non_zero(tab->mat->row[row] + off + tab->n_dead, 2947 tab->n_col - tab->n_dead) == -1; 2948} 2949 2950/* Return the minimal value of the affine expression "f" with denominator 2951 * "denom" in *opt, *opt_denom, assuming the tableau is not empty and 2952 * the expression cannot attain arbitrarily small values. 2953 * If opt_denom is NULL, then *opt is rounded up to the nearest integer. 2954 * The return value reflects the nature of the result (empty, unbounded, 2955 * minimal value returned in *opt). 2956 */ 2957enum isl_lp_result isl_tab_min(struct isl_tab *tab, 2958 isl_int *f, isl_int denom, isl_int *opt, isl_int *opt_denom, 2959 unsigned flags) 2960{ 2961 int r; 2962 enum isl_lp_result res = isl_lp_ok; 2963 struct isl_tab_var *var; 2964 struct isl_tab_undo *snap; 2965 2966 if (!tab) 2967 return isl_lp_error; 2968 2969 if (tab->empty) 2970 return isl_lp_empty; 2971 2972 snap = isl_tab_snap(tab); 2973 r = isl_tab_add_row(tab, f); 2974 if (r < 0) 2975 return isl_lp_error; 2976 var = &tab->con[r]; 2977 for (;;) { 2978 int row, col; 2979 find_pivot(tab, var, var, -1, &row, &col); 2980 if (row == var->index) { 2981 res = isl_lp_unbounded; 2982 break; 2983 } 2984 if (row == -1) 2985 break; 2986 if (isl_tab_pivot(tab, row, col) < 0) 2987 return isl_lp_error; 2988 } 2989 isl_int_mul(tab->mat->row[var->index][0], 2990 tab->mat->row[var->index][0], denom); 2991 if (ISL_FL_ISSET(flags, ISL_TAB_SAVE_DUAL)) { 2992 int i; 2993 2994 isl_vec_free(tab->dual); 2995 tab->dual = isl_vec_alloc(tab->mat->ctx, 1 + tab->n_con); 2996 if (!tab->dual) 2997 return isl_lp_error; 2998 isl_int_set(tab->dual->el[0], tab->mat->row[var->index][0]); 2999 for (i = 0; i < tab->n_con; ++i) { 3000 int pos; 3001 if (tab->con[i].is_row) { 3002 isl_int_set_si(tab->dual->el[1 + i], 0); 3003 continue; 3004 } 3005 pos = 2 + tab->M + tab->con[i].index; 3006 if (tab->con[i].negated) 3007 isl_int_neg(tab->dual->el[1 + i], 3008 tab->mat->row[var->index][pos]); 3009 else 3010 isl_int_set(tab->dual->el[1 + i], 3011 tab->mat->row[var->index][pos]); 3012 } 3013 } 3014 if (opt && res == isl_lp_ok) { 3015 if (opt_denom) { 3016 isl_int_set(*opt, tab->mat->row[var->index][1]); 3017 isl_int_set(*opt_denom, tab->mat->row[var->index][0]); 3018 } else 3019 isl_int_cdiv_q(*opt, tab->mat->row[var->index][1], 3020 tab->mat->row[var->index][0]); 3021 } 3022 if (isl_tab_rollback(tab, snap) < 0) 3023 return isl_lp_error; 3024 return res; 3025} 3026 3027int isl_tab_is_redundant(struct isl_tab *tab, int con) 3028{ 3029 if (!tab) 3030 return -1; 3031 if (tab->con[con].is_zero) 3032 return 0; 3033 if (tab->con[con].is_redundant) 3034 return 1; 3035 return tab->con[con].is_row && tab->con[con].index < tab->n_redundant; 3036} 3037 3038/* Take a snapshot of the tableau that can be restored by s call to 3039 * isl_tab_rollback. 3040 */ 3041struct isl_tab_undo *isl_tab_snap(struct isl_tab *tab) 3042{ 3043 if (!tab) 3044 return NULL; 3045 tab->need_undo = 1; 3046 return tab->top; 3047} 3048 3049/* Undo the operation performed by isl_tab_relax. 3050 */ 3051static int unrelax(struct isl_tab *tab, struct isl_tab_var *var) WARN_UNUSED; 3052static int unrelax(struct isl_tab *tab, struct isl_tab_var *var) 3053{ 3054 unsigned off = 2 + tab->M; 3055 3056 if (!var->is_row && !max_is_manifestly_unbounded(tab, var)) 3057 if (to_row(tab, var, 1) < 0) 3058 return -1; 3059 3060 if (var->is_row) { 3061 isl_int_sub(tab->mat->row[var->index][1], 3062 tab->mat->row[var->index][1], tab->mat->row[var->index][0]); 3063 if (var->is_nonneg) { 3064 int sgn = restore_row(tab, var); 3065 isl_assert(tab->mat->ctx, sgn >= 0, return -1); 3066 } 3067 } else { 3068 int i; 3069 3070 for (i = 0; i < tab->n_row; ++i) { 3071 if (isl_int_is_zero(tab->mat->row[i][off + var->index])) 3072 continue; 3073 isl_int_add(tab->mat->row[i][1], tab->mat->row[i][1], 3074 tab->mat->row[i][off + var->index]); 3075 } 3076 3077 } 3078 3079 return 0; 3080} 3081 3082/* Undo the operation performed by isl_tab_unrestrict. 3083 * 3084 * In particular, mark the variable as being non-negative and make 3085 * sure the sample value respects this constraint. 3086 */ 3087static int ununrestrict(struct isl_tab *tab, struct isl_tab_var *var) 3088{ 3089 var->is_nonneg = 1; 3090 3091 if (var->is_row && restore_row(tab, var) < -1) 3092 return -1; 3093 3094 return 0; 3095} 3096 3097static int perform_undo_var(struct isl_tab *tab, struct isl_tab_undo *undo) WARN_UNUSED; 3098static int perform_undo_var(struct isl_tab *tab, struct isl_tab_undo *undo) 3099{ 3100 struct isl_tab_var *var = var_from_index(tab, undo->u.var_index); 3101 switch (undo->type) { 3102 case isl_tab_undo_nonneg: 3103 var->is_nonneg = 0; 3104 break; 3105 case isl_tab_undo_redundant: 3106 var->is_redundant = 0; 3107 tab->n_redundant--; 3108 restore_row(tab, isl_tab_var_from_row(tab, tab->n_redundant)); 3109 break; 3110 case isl_tab_undo_freeze: 3111 var->frozen = 0; 3112 break; 3113 case isl_tab_undo_zero: 3114 var->is_zero = 0; 3115 if (!var->is_row) 3116 tab->n_dead--; 3117 break; 3118 case isl_tab_undo_allocate: 3119 if (undo->u.var_index >= 0) { 3120 isl_assert(tab->mat->ctx, !var->is_row, return -1); 3121 drop_col(tab, var->index); 3122 break; 3123 } 3124 if (!var->is_row) { 3125 if (!max_is_manifestly_unbounded(tab, var)) { 3126 if (to_row(tab, var, 1) < 0) 3127 return -1; 3128 } else if (!min_is_manifestly_unbounded(tab, var)) { 3129 if (to_row(tab, var, -1) < 0) 3130 return -1; 3131 } else 3132 if (to_row(tab, var, 0) < 0) 3133 return -1; 3134 } 3135 drop_row(tab, var->index); 3136 break; 3137 case isl_tab_undo_relax: 3138 return unrelax(tab, var); 3139 case isl_tab_undo_unrestrict: 3140 return ununrestrict(tab, var); 3141 default: 3142 isl_die(tab->mat->ctx, isl_error_internal, 3143 "perform_undo_var called on invalid undo record", 3144 return -1); 3145 } 3146 3147 return 0; 3148} 3149 3150/* Restore the tableau to the state where the basic variables 3151 * are those in "col_var". 3152 * We first construct a list of variables that are currently in 3153 * the basis, but shouldn't. Then we iterate over all variables 3154 * that should be in the basis and for each one that is currently 3155 * not in the basis, we exchange it with one of the elements of the 3156 * list constructed before. 3157 * We can always find an appropriate variable to pivot with because 3158 * the current basis is mapped to the old basis by a non-singular 3159 * matrix and so we can never end up with a zero row. 3160 */ 3161static int restore_basis(struct isl_tab *tab, int *col_var) 3162{ 3163 int i, j; 3164 int n_extra = 0; 3165 int *extra = NULL; /* current columns that contain bad stuff */ 3166 unsigned off = 2 + tab->M; 3167 3168 extra = isl_alloc_array(tab->mat->ctx, int, tab->n_col); 3169 if (tab->n_col && !extra) 3170 goto error; 3171 for (i = 0; i < tab->n_col; ++i) { 3172 for (j = 0; j < tab->n_col; ++j) 3173 if (tab->col_var[i] == col_var[j]) 3174 break; 3175 if (j < tab->n_col) 3176 continue; 3177 extra[n_extra++] = i; 3178 } 3179 for (i = 0; i < tab->n_col && n_extra > 0; ++i) { 3180 struct isl_tab_var *var; 3181 int row; 3182 3183 for (j = 0; j < tab->n_col; ++j) 3184 if (col_var[i] == tab->col_var[j]) 3185 break; 3186 if (j < tab->n_col) 3187 continue; 3188 var = var_from_index(tab, col_var[i]); 3189 row = var->index; 3190 for (j = 0; j < n_extra; ++j) 3191 if (!isl_int_is_zero(tab->mat->row[row][off+extra[j]])) 3192 break; 3193 isl_assert(tab->mat->ctx, j < n_extra, goto error); 3194 if (isl_tab_pivot(tab, row, extra[j]) < 0) 3195 goto error; 3196 extra[j] = extra[--n_extra]; 3197 } 3198 3199 free(extra); 3200 return 0; 3201error: 3202 free(extra); 3203 return -1; 3204} 3205 3206/* Remove all samples with index n or greater, i.e., those samples 3207 * that were added since we saved this number of samples in 3208 * isl_tab_save_samples. 3209 */ 3210static void drop_samples_since(struct isl_tab *tab, int n) 3211{ 3212 int i; 3213 3214 for (i = tab->n_sample - 1; i >= 0 && tab->n_sample > n; --i) { 3215 if (tab->sample_index[i] < n) 3216 continue; 3217 3218 if (i != tab->n_sample - 1) { 3219 int t = tab->sample_index[tab->n_sample-1]; 3220 tab->sample_index[tab->n_sample-1] = tab->sample_index[i]; 3221 tab->sample_index[i] = t; 3222 isl_mat_swap_rows(tab->samples, tab->n_sample-1, i); 3223 } 3224 tab->n_sample--; 3225 } 3226} 3227 3228static int perform_undo(struct isl_tab *tab, struct isl_tab_undo *undo) WARN_UNUSED; 3229static int perform_undo(struct isl_tab *tab, struct isl_tab_undo *undo) 3230{ 3231 switch (undo->type) { 3232 case isl_tab_undo_empty: 3233 tab->empty = 0; 3234 break; 3235 case isl_tab_undo_nonneg: 3236 case isl_tab_undo_redundant: 3237 case isl_tab_undo_freeze: 3238 case isl_tab_undo_zero: 3239 case isl_tab_undo_allocate: 3240 case isl_tab_undo_relax: 3241 case isl_tab_undo_unrestrict: 3242 return perform_undo_var(tab, undo); 3243 case isl_tab_undo_bmap_eq: 3244 return isl_basic_map_free_equality(tab->bmap, 1); 3245 case isl_tab_undo_bmap_ineq: 3246 return isl_basic_map_free_inequality(tab->bmap, 1); 3247 case isl_tab_undo_bmap_div: 3248 if (isl_basic_map_free_div(tab->bmap, 1) < 0) 3249 return -1; 3250 if (tab->samples) 3251 tab->samples->n_col--; 3252 break; 3253 case isl_tab_undo_saved_basis: 3254 if (restore_basis(tab, undo->u.col_var) < 0) 3255 return -1; 3256 break; 3257 case isl_tab_undo_drop_sample: 3258 tab->n_outside--; 3259 break; 3260 case isl_tab_undo_saved_samples: 3261 drop_samples_since(tab, undo->u.n); 3262 break; 3263 case isl_tab_undo_callback: 3264 return undo->u.callback->run(undo->u.callback); 3265 default: 3266 isl_assert(tab->mat->ctx, 0, return -1); 3267 } 3268 return 0; 3269} 3270 3271/* Return the tableau to the state it was in when the snapshot "snap" 3272 * was taken. 3273 */ 3274int isl_tab_rollback(struct isl_tab *tab, struct isl_tab_undo *snap) 3275{ 3276 struct isl_tab_undo *undo, *next; 3277 3278 if (!tab) 3279 return -1; 3280 3281 tab->in_undo = 1; 3282 for (undo = tab->top; undo && undo != &tab->bottom; undo = next) { 3283 next = undo->next; 3284 if (undo == snap) 3285 break; 3286 if (perform_undo(tab, undo) < 0) { 3287 tab->top = undo; 3288 free_undo(tab); 3289 tab->in_undo = 0; 3290 return -1; 3291 } 3292 free_undo_record(undo); 3293 } 3294 tab->in_undo = 0; 3295 tab->top = undo; 3296 if (!undo) 3297 return -1; 3298 return 0; 3299} 3300 3301/* The given row "row" represents an inequality violated by all 3302 * points in the tableau. Check for some special cases of such 3303 * separating constraints. 3304 * In particular, if the row has been reduced to the constant -1, 3305 * then we know the inequality is adjacent (but opposite) to 3306 * an equality in the tableau. 3307 * If the row has been reduced to r = c*(-1 -r'), with r' an inequality 3308 * of the tableau and c a positive constant, then the inequality 3309 * is adjacent (but opposite) to the inequality r'. 3310 */ 3311static enum isl_ineq_type separation_type(struct isl_tab *tab, unsigned row) 3312{ 3313 int pos; 3314 unsigned off = 2 + tab->M; 3315 3316 if (tab->rational) 3317 return isl_ineq_separate; 3318 3319 if (!isl_int_is_one(tab->mat->row[row][0])) 3320 return isl_ineq_separate; 3321 3322 pos = isl_seq_first_non_zero(tab->mat->row[row] + off + tab->n_dead, 3323 tab->n_col - tab->n_dead); 3324 if (pos == -1) { 3325 if (isl_int_is_negone(tab->mat->row[row][1])) 3326 return isl_ineq_adj_eq; 3327 else 3328 return isl_ineq_separate; 3329 } 3330 3331 if (!isl_int_eq(tab->mat->row[row][1], 3332 tab->mat->row[row][off + tab->n_dead + pos])) 3333 return isl_ineq_separate; 3334 3335 pos = isl_seq_first_non_zero( 3336 tab->mat->row[row] + off + tab->n_dead + pos + 1, 3337 tab->n_col - tab->n_dead - pos - 1); 3338 3339 return pos == -1 ? isl_ineq_adj_ineq : isl_ineq_separate; 3340} 3341 3342/* Check the effect of inequality "ineq" on the tableau "tab". 3343 * The result may be 3344 * isl_ineq_redundant: satisfied by all points in the tableau 3345 * isl_ineq_separate: satisfied by no point in the tableau 3346 * isl_ineq_cut: satisfied by some by not all points 3347 * isl_ineq_adj_eq: adjacent to an equality 3348 * isl_ineq_adj_ineq: adjacent to an inequality. 3349 */ 3350enum isl_ineq_type isl_tab_ineq_type(struct isl_tab *tab, isl_int *ineq) 3351{ 3352 enum isl_ineq_type type = isl_ineq_error; 3353 struct isl_tab_undo *snap = NULL; 3354 int con; 3355 int row; 3356 3357 if (!tab) 3358 return isl_ineq_error; 3359 3360 if (isl_tab_extend_cons(tab, 1) < 0) 3361 return isl_ineq_error; 3362 3363 snap = isl_tab_snap(tab); 3364 3365 con = isl_tab_add_row(tab, ineq); 3366 if (con < 0) 3367 goto error; 3368 3369 row = tab->con[con].index; 3370 if (isl_tab_row_is_redundant(tab, row)) 3371 type = isl_ineq_redundant; 3372 else if (isl_int_is_neg(tab->mat->row[row][1]) && 3373 (tab->rational || 3374 isl_int_abs_ge(tab->mat->row[row][1], 3375 tab->mat->row[row][0]))) { 3376 int nonneg = at_least_zero(tab, &tab->con[con]); 3377 if (nonneg < 0) 3378 goto error; 3379 if (nonneg) 3380 type = isl_ineq_cut; 3381 else 3382 type = separation_type(tab, row); 3383 } else { 3384 int red = con_is_redundant(tab, &tab->con[con]); 3385 if (red < 0) 3386 goto error; 3387 if (!red) 3388 type = isl_ineq_cut; 3389 else 3390 type = isl_ineq_redundant; 3391 } 3392 3393 if (isl_tab_rollback(tab, snap)) 3394 return isl_ineq_error; 3395 return type; 3396error: 3397 return isl_ineq_error; 3398} 3399 3400int isl_tab_track_bmap(struct isl_tab *tab, __isl_take isl_basic_map *bmap) 3401{ 3402 bmap = isl_basic_map_cow(bmap); 3403 if (!tab || !bmap) 3404 goto error; 3405 3406 if (tab->empty) { 3407 bmap = isl_basic_map_set_to_empty(bmap); 3408 if (!bmap) 3409 goto error; 3410 tab->bmap = bmap; 3411 return 0; 3412 } 3413 3414 isl_assert(tab->mat->ctx, tab->n_eq == bmap->n_eq, goto error); 3415 isl_assert(tab->mat->ctx, 3416 tab->n_con == bmap->n_eq + bmap->n_ineq, goto error); 3417 3418 tab->bmap = bmap; 3419 3420 return 0; 3421error: 3422 isl_basic_map_free(bmap); 3423 return -1; 3424} 3425 3426int isl_tab_track_bset(struct isl_tab *tab, __isl_take isl_basic_set *bset) 3427{ 3428 return isl_tab_track_bmap(tab, (isl_basic_map *)bset); 3429} 3430 3431__isl_keep isl_basic_set *isl_tab_peek_bset(struct isl_tab *tab) 3432{ 3433 if (!tab) 3434 return NULL; 3435 3436 return (isl_basic_set *)tab->bmap; 3437} 3438 3439static void isl_tab_print_internal(__isl_keep struct isl_tab *tab, 3440 FILE *out, int indent) 3441{ 3442 unsigned r, c; 3443 int i; 3444 3445 if (!tab) { 3446 fprintf(out, "%*snull tab\n", indent, ""); 3447 return; 3448 } 3449 fprintf(out, "%*sn_redundant: %d, n_dead: %d", indent, "", 3450 tab->n_redundant, tab->n_dead); 3451 if (tab->rational) 3452 fprintf(out, ", rational"); 3453 if (tab->empty) 3454 fprintf(out, ", empty"); 3455 fprintf(out, "\n"); 3456 fprintf(out, "%*s[", indent, ""); 3457 for (i = 0; i < tab->n_var; ++i) { 3458 if (i) 3459 fprintf(out, (i == tab->n_param || 3460 i == tab->n_var - tab->n_div) ? "; " 3461 : ", "); 3462 fprintf(out, "%c%d%s", tab->var[i].is_row ? 'r' : 'c', 3463 tab->var[i].index, 3464 tab->var[i].is_zero ? " [=0]" : 3465 tab->var[i].is_redundant ? " [R]" : ""); 3466 } 3467 fprintf(out, "]\n"); 3468 fprintf(out, "%*s[", indent, ""); 3469 for (i = 0; i < tab->n_con; ++i) { 3470 if (i) 3471 fprintf(out, ", "); 3472 fprintf(out, "%c%d%s", tab->con[i].is_row ? 'r' : 'c', 3473 tab->con[i].index, 3474 tab->con[i].is_zero ? " [=0]" : 3475 tab->con[i].is_redundant ? " [R]" : ""); 3476 } 3477 fprintf(out, "]\n"); 3478 fprintf(out, "%*s[", indent, ""); 3479 for (i = 0; i < tab->n_row; ++i) { 3480 const char *sign = ""; 3481 if (i) 3482 fprintf(out, ", "); 3483 if (tab->row_sign) { 3484 if (tab->row_sign[i] == isl_tab_row_unknown) 3485 sign = "?"; 3486 else if (tab->row_sign[i] == isl_tab_row_neg) 3487 sign = "-"; 3488 else if (tab->row_sign[i] == isl_tab_row_pos) 3489 sign = "+"; 3490 else 3491 sign = "+-"; 3492 } 3493 fprintf(out, "r%d: %d%s%s", i, tab->row_var[i], 3494 isl_tab_var_from_row(tab, i)->is_nonneg ? " [>=0]" : "", sign); 3495 } 3496 fprintf(out, "]\n"); 3497 fprintf(out, "%*s[", indent, ""); 3498 for (i = 0; i < tab->n_col; ++i) { 3499 if (i) 3500 fprintf(out, ", "); 3501 fprintf(out, "c%d: %d%s", i, tab->col_var[i], 3502 var_from_col(tab, i)->is_nonneg ? " [>=0]" : ""); 3503 } 3504 fprintf(out, "]\n"); 3505 r = tab->mat->n_row; 3506 tab->mat->n_row = tab->n_row; 3507 c = tab->mat->n_col; 3508 tab->mat->n_col = 2 + tab->M + tab->n_col; 3509 isl_mat_print_internal(tab->mat, out, indent); 3510 tab->mat->n_row = r; 3511 tab->mat->n_col = c; 3512 if (tab->bmap) 3513 isl_basic_map_print_internal(tab->bmap, out, indent); 3514} 3515 3516void isl_tab_dump(__isl_keep struct isl_tab *tab) 3517{ 3518 isl_tab_print_internal(tab, stderr, 0); 3519} 3520