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