Coverage Report

Created: 2017-03-28 09:59

/Users/buildslave/jenkins/sharedspace/clang-stage2-coverage-R@2/llvm/tools/polly/lib/External/isl/isl_tab_pip.c
Line
Count
Source (jump to first uncovered line)
1
/*
2
 * Copyright 2008-2009 Katholieke Universiteit Leuven
3
 * Copyright 2010      INRIA Saclay
4
 * Copyright 2016      Sven Verdoolaege
5
 *
6
 * Use of this software is governed by the MIT license
7
 *
8
 * Written by Sven Verdoolaege, K.U.Leuven, Departement
9
 * Computerwetenschappen, Celestijnenlaan 200A, B-3001 Leuven, Belgium
10
 * and INRIA Saclay - Ile-de-France, Parc Club Orsay Universite,
11
 * ZAC des vignes, 4 rue Jacques Monod, 91893 Orsay, France 
12
 */
13
14
#include <isl_ctx_private.h>
15
#include "isl_map_private.h"
16
#include <isl_seq.h>
17
#include "isl_tab.h"
18
#include "isl_sample.h"
19
#include <isl_mat_private.h>
20
#include <isl_vec_private.h>
21
#include <isl_aff_private.h>
22
#include <isl_constraint_private.h>
23
#include <isl_options_private.h>
24
#include <isl_config.h>
25
26
#include <bset_to_bmap.c>
27
28
/*
29
 * The implementation of parametric integer linear programming in this file
30
 * was inspired by the paper "Parametric Integer Programming" and the
31
 * report "Solving systems of affine (in)equalities" by Paul Feautrier
32
 * (and others).
33
 *
34
 * The strategy used for obtaining a feasible solution is different
35
 * from the one used in isl_tab.c.  In particular, in isl_tab.c,
36
 * upon finding a constraint that is not yet satisfied, we pivot
37
 * in a row that increases the constant term of the row holding the
38
 * constraint, making sure the sample solution remains feasible
39
 * for all the constraints it already satisfied.
40
 * Here, we always pivot in the row holding the constraint,
41
 * choosing a column that induces the lexicographically smallest
42
 * increment to the sample solution.
43
 *
44
 * By starting out from a sample value that is lexicographically
45
 * smaller than any integer point in the problem space, the first
46
 * feasible integer sample point we find will also be the lexicographically
47
 * smallest.  If all variables can be assumed to be non-negative,
48
 * then the initial sample value may be chosen equal to zero.
49
 * However, we will not make this assumption.  Instead, we apply
50
 * the "big parameter" trick.  Any variable x is then not directly
51
 * used in the tableau, but instead it is represented by another
52
 * variable x' = M + x, where M is an arbitrarily large (positive)
53
 * value.  x' is therefore always non-negative, whatever the value of x.
54
 * Taking as initial sample value x' = 0 corresponds to x = -M,
55
 * which is always smaller than any possible value of x.
56
 *
57
 * The big parameter trick is used in the main tableau and
58
 * also in the context tableau if isl_context_lex is used.
59
 * In this case, each tableaus has its own big parameter.
60
 * Before doing any real work, we check if all the parameters
61
 * happen to be non-negative.  If so, we drop the column corresponding
62
 * to M from the initial context tableau.
63
 * If isl_context_gbr is used, then the big parameter trick is only
64
 * used in the main tableau.
65
 */
66
67
struct isl_context;
68
struct isl_context_op {
69
  /* detect nonnegative parameters in context and mark them in tab */
70
  struct isl_tab *(*detect_nonnegative_parameters)(
71
      struct isl_context *context, struct isl_tab *tab);
72
  /* return temporary reference to basic set representation of context */
73
  struct isl_basic_set *(*peek_basic_set)(struct isl_context *context);
74
  /* return temporary reference to tableau representation of context */
75
  struct isl_tab *(*peek_tab)(struct isl_context *context);
76
  /* add equality; check is 1 if eq may not be valid;
77
   * update is 1 if we may want to call ineq_sign on context later.
78
   */
79
  void (*add_eq)(struct isl_context *context, isl_int *eq,
80
      int check, int update);
81
  /* add inequality; check is 1 if ineq may not be valid;
82
   * update is 1 if we may want to call ineq_sign on context later.
83
   */
84
  void (*add_ineq)(struct isl_context *context, isl_int *ineq,
85
      int check, int update);
86
  /* check sign of ineq based on previous information.
87
   * strict is 1 if saturation should be treated as a positive sign.
88
   */
89
  enum isl_tab_row_sign (*ineq_sign)(struct isl_context *context,
90
      isl_int *ineq, int strict);
91
  /* check if inequality maintains feasibility */
92
  int (*test_ineq)(struct isl_context *context, isl_int *ineq);
93
  /* return index of a div that corresponds to "div" */
94
  int (*get_div)(struct isl_context *context, struct isl_tab *tab,
95
      struct isl_vec *div);
96
  /* insert div "div" to context at "pos" and return non-negativity */
97
  isl_bool (*insert_div)(struct isl_context *context, int pos,
98
    __isl_keep isl_vec *div);
99
  int (*detect_equalities)(struct isl_context *context,
100
      struct isl_tab *tab);
101
  /* return row index of "best" split */
102
  int (*best_split)(struct isl_context *context, struct isl_tab *tab);
103
  /* check if context has already been determined to be empty */
104
  int (*is_empty)(struct isl_context *context);
105
  /* check if context is still usable */
106
  int (*is_ok)(struct isl_context *context);
107
  /* save a copy/snapshot of context */
108
  void *(*save)(struct isl_context *context);
109
  /* restore saved context */
110
  void (*restore)(struct isl_context *context, void *);
111
  /* discard saved context */
112
  void (*discard)(void *);
113
  /* invalidate context */
114
  void (*invalidate)(struct isl_context *context);
115
  /* free context */
116
  __isl_null struct isl_context *(*free)(struct isl_context *context);
117
};
118
119
/* Shared parts of context representation.
120
 *
121
 * "n_unknown" is the number of final unknown integer divisions
122
 * in the input domain.
123
 */
124
struct isl_context {
125
  struct isl_context_op *op;
126
  int n_unknown;
127
};
128
129
struct isl_context_lex {
130
  struct isl_context context;
131
  struct isl_tab *tab;
132
};
133
134
/* A stack (linked list) of solutions of subtrees of the search space.
135
 *
136
 * "ma" describes the solution as a function of "dom".
137
 * In particular, the domain space of "ma" is equal to the space of "dom".
138
 *
139
 * If "ma" is NULL, then there is no solution on "dom".
140
 */
141
struct isl_partial_sol {
142
  int level;
143
  struct isl_basic_set *dom;
144
  isl_multi_aff *ma;
145
146
  struct isl_partial_sol *next;
147
};
148
149
struct isl_sol;
150
struct isl_sol_callback {
151
  struct isl_tab_callback callback;
152
  struct isl_sol *sol;
153
};
154
155
/* isl_sol is an interface for constructing a solution to
156
 * a parametric integer linear programming problem.
157
 * Every time the algorithm reaches a state where a solution
158
 * can be read off from the tableau, the function "add" is called
159
 * on the isl_sol passed to find_solutions_main.  In a state where
160
 * the tableau is empty, "add_empty" is called instead.
161
 * "free" is called to free the implementation specific fields, if any.
162
 *
163
 * "error" is set if some error has occurred.  This flag invalidates
164
 * the remainder of the data structure.
165
 * If "rational" is set, then a rational optimization is being performed.
166
 * "level" is the current level in the tree with nodes for each
167
 * split in the context.
168
 * If "max" is set, then a maximization problem is being solved, rather than
169
 * a minimization problem, which means that the variables in the
170
 * tableau have value "M - x" rather than "M + x".
171
 * "n_out" is the number of output dimensions in the input.
172
 * "space" is the space in which the solution (and also the input) lives.
173
 *
174
 * The context tableau is owned by isl_sol and is updated incrementally.
175
 *
176
 * There are currently three implementations of this interface,
177
 * isl_sol_map, which simply collects the solutions in an isl_map
178
 * and (optionally) the parts of the context where there is no solution
179
 * in an isl_set,
180
 * isl_sol_pma, which collects an isl_pw_multi_aff instead, and
181
 * isl_sol_for, which calls a user-defined function for each part of
182
 * the solution.
183
 */
184
struct isl_sol {
185
  int error;
186
  int rational;
187
  int level;
188
  int max;
189
  int n_out;
190
  isl_space *space;
191
  struct isl_context *context;
192
  struct isl_partial_sol *partial;
193
  void (*add)(struct isl_sol *sol,
194
    __isl_take isl_basic_set *dom, __isl_take isl_multi_aff *ma);
195
  void (*add_empty)(struct isl_sol *sol, struct isl_basic_set *bset);
196
  void (*free)(struct isl_sol *sol);
197
  struct isl_sol_callback dec_level;
198
};
199
200
static void sol_free(struct isl_sol *sol)
201
3.22k
{
202
3.22k
  struct isl_partial_sol *partial, *next;
203
3.22k
  if (!sol)
204
0
    return;
205
3.22k
  
for (partial = sol->partial; 3.22k
partial3.22k
;
partial = next0
)
{0
206
0
    next = partial->next;
207
0
    isl_basic_set_free(partial->dom);
208
0
    isl_multi_aff_free(partial->ma);
209
0
    free(partial);
210
0
  }
211
3.22k
  isl_space_free(sol->space);
212
3.22k
  if (sol->context)
213
3.22k
    sol->context->op->free(sol->context);
214
3.22k
  sol->free(sol);
215
3.22k
  free(sol);
216
3.22k
}
217
218
/* Push a partial solution represented by a domain and function "ma"
219
 * onto the stack of partial solutions.
220
 * If "ma" is NULL, then "dom" represents a part of the domain
221
 * with no solution.
222
 */
223
static void sol_push_sol(struct isl_sol *sol,
224
  __isl_take isl_basic_set *dom, __isl_take isl_multi_aff *ma)
225
3.87k
{
226
3.87k
  struct isl_partial_sol *partial;
227
3.87k
228
3.87k
  if (
sol->error || 3.87k
!dom3.87k
)
229
0
    goto error;
230
3.87k
231
3.87k
  
partial = 3.87k
isl_alloc_type3.87k
(dom->ctx, struct isl_partial_sol);
232
3.87k
  if (!partial)
233
0
    goto error;
234
3.87k
235
3.87k
  partial->level = sol->level;
236
3.87k
  partial->dom = dom;
237
3.87k
  partial->ma = ma;
238
3.87k
  partial->next = sol->partial;
239
3.87k
240
3.87k
  sol->partial = partial;
241
3.87k
242
3.87k
  return;
243
0
error:
244
0
  isl_basic_set_free(dom);
245
0
  isl_multi_aff_free(ma);
246
0
  sol->error = 1;
247
0
}
248
249
/* Check that the final columns of "M", starting at "first", are zero.
250
 */
251
static isl_stat check_final_columns_are_zero(__isl_keep isl_mat *M,
252
  unsigned first)
253
2.95k
{
254
2.95k
  int i;
255
2.95k
  unsigned rows, cols, n;
256
2.95k
257
2.95k
  if (!M)
258
0
    return isl_stat_error;
259
2.95k
  rows = isl_mat_rows(M);
260
2.95k
  cols = isl_mat_cols(M);
261
2.95k
  n = cols - first;
262
12.1k
  for (i = 0; 
i < rows12.1k
;
++i9.19k
)
263
9.19k
    
if (9.19k
isl_seq_first_non_zero(M->row[i] + first, n) != -19.19k
)
264
0
      isl_die(isl_mat_get_ctx(M), isl_error_internal,
265
2.95k
        "final columns should be zero",
266
2.95k
        return isl_stat_error);
267
2.95k
  return isl_stat_ok;
268
2.95k
}
269
270
/* Set the affine expressions in "ma" according to the rows in "M", which
271
 * are defined over the local space "ls".
272
 * The matrix "M" may have extra (zero) columns beyond the number
273
 * of variables in "ls".
274
 */
275
static __isl_give isl_multi_aff *set_from_affine_matrix(
276
  __isl_take isl_multi_aff *ma, __isl_take isl_local_space *ls,
277
  __isl_take isl_mat *M)
278
2.95k
{
279
2.95k
  int i, dim;
280
2.95k
  isl_aff *aff;
281
2.95k
282
2.95k
  if (
!ma || 2.95k
!ls2.95k
||
!M2.95k
)
283
0
    goto error;
284
2.95k
285
2.95k
  dim = isl_local_space_dim(ls, isl_dim_all);
286
2.95k
  if (check_final_columns_are_zero(M, 1 + dim) < 0)
287
0
    goto error;
288
9.19k
  
for (i = 1; 2.95k
i < M->n_row9.19k
;
++i6.23k
)
{6.23k
289
6.23k
    aff = isl_aff_alloc(isl_local_space_copy(ls));
290
6.23k
    if (
aff6.23k
)
{6.23k
291
6.23k
      isl_int_set(aff->v->el[0], M->row[0][0]);
292
6.23k
      isl_seq_cpy(aff->v->el + 1, M->row[i], 1 + dim);
293
6.23k
    }
294
6.23k
    aff = isl_aff_normalize(aff);
295
6.23k
    ma = isl_multi_aff_set_aff(ma, i - 1, aff);
296
6.23k
  }
297
2.95k
  isl_local_space_free(ls);
298
2.95k
  isl_mat_free(M);
299
2.95k
300
2.95k
  return ma;
301
0
error:
302
0
  isl_local_space_free(ls);
303
0
  isl_mat_free(M);
304
0
  isl_multi_aff_free(ma);
305
0
  return NULL;
306
2.95k
}
307
308
/* Push a partial solution represented by a domain and mapping M
309
 * onto the stack of partial solutions.
310
 *
311
 * The affine matrix "M" maps the dimensions of the context
312
 * to the output variables.  Convert it into an isl_multi_aff and
313
 * then call sol_push_sol.
314
 *
315
 * Note that the description of the initial context may have involved
316
 * existentially quantified variables, in which case they also appear
317
 * in "dom".  These need to be removed before creating the affine
318
 * expression because an affine expression cannot be defined in terms
319
 * of existentially quantified variables without a known representation.
320
 * Since newly added integer divisions are inserted before these
321
 * existentially quantified variables, they are still in the final
322
 * positions and the corresponding final columns of "M" are zero
323
 * because align_context_divs adds the existentially quantified
324
 * variables of the context to the main tableau without any constraints and
325
 * any equality constraints that are added later on can only serve
326
 * to eliminate these existentially quantified variables.
327
 */
328
static void sol_push_sol_mat(struct isl_sol *sol,
329
  __isl_take isl_basic_set *dom, __isl_take isl_mat *M)
330
2.95k
{
331
2.95k
  isl_local_space *ls;
332
2.95k
  isl_multi_aff *ma;
333
2.95k
  int n_div, n_known;
334
2.95k
335
2.95k
  n_div = isl_basic_set_dim(dom, isl_dim_div);
336
2.95k
  n_known = n_div - sol->context->n_unknown;
337
2.95k
338
2.95k
  ma = isl_multi_aff_alloc(isl_space_copy(sol->space));
339
2.95k
  ls = isl_basic_set_get_local_space(dom);
340
2.95k
  ls = isl_local_space_drop_dims(ls, isl_dim_div,
341
2.95k
          n_known, n_div - n_known);
342
2.95k
  ma = set_from_affine_matrix(ma, ls, M);
343
2.95k
344
2.95k
  if (!ma)
345
0
    dom = isl_basic_set_free(dom);
346
2.95k
  sol_push_sol(sol, dom, ma);
347
2.95k
}
348
349
/* Pop one partial solution from the partial solution stack and
350
 * pass it on to sol->add or sol->add_empty.
351
 */
352
static void sol_pop_one(struct isl_sol *sol)
353
3.79k
{
354
3.79k
  struct isl_partial_sol *partial;
355
3.79k
356
3.79k
  partial = sol->partial;
357
3.79k
  sol->partial = partial->next;
358
3.79k
359
3.79k
  if (partial->ma)
360
2.87k
    sol->add(sol, partial->dom, partial->ma);
361
3.79k
  else
362
919
    sol->add_empty(sol, partial->dom);
363
3.79k
  free(partial);
364
3.79k
}
365
366
/* Return a fresh copy of the domain represented by the context tableau.
367
 */
368
static struct isl_basic_set *sol_domain(struct isl_sol *sol)
369
3.96k
{
370
3.96k
  struct isl_basic_set *bset;
371
3.96k
372
3.96k
  if (sol->error)
373
0
    return NULL;
374
3.96k
375
3.96k
  bset = isl_basic_set_dup(sol->context->op->peek_basic_set(sol->context));
376
3.96k
  bset = isl_basic_set_update_from_tab(bset,
377
3.96k
      sol->context->op->peek_tab(sol->context));
378
3.96k
379
3.96k
  return bset;
380
3.96k
}
381
382
/* Check whether two partial solutions have the same affine expressions.
383
 */
384
static isl_bool same_solution(struct isl_partial_sol *s1,
385
  struct isl_partial_sol *s2)
386
772
{
387
772
  if (!s1->ma != !s2->ma)
388
676
    return isl_bool_false;
389
96
  
if (96
!s1->ma96
)
390
0
    return isl_bool_true;
391
96
392
96
  return isl_multi_aff_plain_is_equal(s1->ma, s2->ma);
393
96
}
394
395
/* Swap the initial two partial solutions in "sol".
396
 *
397
 * That is, go from
398
 *
399
 *  sol->partial = p1; p1->next = p2; p2->next = p3
400
 *
401
 * to
402
 *
403
 *  sol->partial = p2; p2->next = p1; p1->next = p3
404
 */
405
static void swap_initial(struct isl_sol *sol)
406
36
{
407
36
  struct isl_partial_sol *partial;
408
36
409
36
  partial = sol->partial;
410
36
  sol->partial = partial->next;
411
36
  partial->next = partial->next->next;
412
36
  sol->partial->next = partial;
413
36
}
414
415
/* Combine the initial two partial solution of "sol" into
416
 * a partial solution with the current context domain of "sol" and
417
 * the function description of the second partial solution in the list.
418
 * The level of the new partial solution is set to the current level.
419
 *
420
 * That is, the first two partial solutions (D1,M1) and (D2,M2) are
421
 * replaced by (D,M2), where D is the domain of "sol", which is assumed
422
 * to be the union of D1 and D2, while M1 is assumed to be equal to M2
423
 * (at least on D1).
424
 */
425
static isl_stat combine_initial_into_second(struct isl_sol *sol)
426
82
{
427
82
  struct isl_partial_sol *partial;
428
82
  isl_basic_set *bset;
429
82
430
82
  partial = sol->partial;
431
82
432
82
  bset = sol_domain(sol);
433
82
  isl_basic_set_free(partial->next->dom);
434
82
  partial->next->dom = bset;
435
82
  partial->next->level = sol->level;
436
82
437
82
  if (!bset)
438
0
    return isl_stat_error;
439
82
440
82
  sol->partial = partial->next;
441
82
  isl_basic_set_free(partial->dom);
442
82
  isl_multi_aff_free(partial->ma);
443
82
  free(partial);
444
82
445
82
  return isl_stat_ok;
446
82
}
447
448
/* Are "ma1" and "ma2" equal to each other on "dom"?
449
 *
450
 * Combine "ma1" and "ma2" with "dom" and check if the results are the same.
451
 * "dom" may have existentially quantified variables.  Eliminate them first
452
 * as otherwise they would have to be eliminated twice, in a more complicated
453
 * context.
454
 */
455
static isl_bool equal_on_domain(__isl_keep isl_multi_aff *ma1,
456
  __isl_keep isl_multi_aff *ma2, __isl_keep isl_basic_set *dom)
457
137
{
458
137
  isl_set *set;
459
137
  isl_pw_multi_aff *pma1, *pma2;
460
137
  isl_bool equal;
461
137
462
137
  set = isl_basic_set_compute_divs(isl_basic_set_copy(dom));
463
137
  pma1 = isl_pw_multi_aff_alloc(isl_set_copy(set),
464
137
          isl_multi_aff_copy(ma1));
465
137
  pma2 = isl_pw_multi_aff_alloc(set, isl_multi_aff_copy(ma2));
466
137
  equal = isl_pw_multi_aff_is_equal(pma1, pma2);
467
137
  isl_pw_multi_aff_free(pma1);
468
137
  isl_pw_multi_aff_free(pma2);
469
137
470
137
  return equal;
471
137
}
472
473
/* The initial two partial solutions of "sol" are known to be at
474
 * the same level.
475
 * If they represent the same solution (on different parts of the domain),
476
 * then combine them into a single solution at the current level.
477
 * Otherwise, pop them both.
478
 *
479
 * Even if the two partial solution are not obviously the same,
480
 * one may still be a simplification of the other over its own domain.
481
 * Also check if the two sets of affine functions are equal when
482
 * restricted to one of the domains.  If so, combine the two
483
 * using the set of affine functions on the other domain.
484
 * That is, for two partial solutions (D1,M1) and (D2,M2),
485
 * if M1 = M2 on D1, then the pair of partial solutions can
486
 * be replaced by (D1+D2,M2) and similarly when M1 = M2 on D2.
487
 */
488
static isl_stat combine_initial_if_equal(struct isl_sol *sol)
489
772
{
490
772
  struct isl_partial_sol *partial;
491
772
  isl_bool same;
492
772
493
772
  partial = sol->partial;
494
772
495
772
  same = same_solution(partial, partial->next);
496
772
  if (same < 0)
497
0
    return isl_stat_error;
498
772
  
if (772
same772
)
499
9
    return combine_initial_into_second(sol);
500
763
  
if (763
partial->ma && 763
partial->next->ma97
)
{87
501
87
    same = equal_on_domain(partial->ma, partial->next->ma,
502
87
          partial->dom);
503
87
    if (same < 0)
504
0
      return isl_stat_error;
505
87
    
if (87
same87
)
506
37
      return combine_initial_into_second(sol);
507
50
    same = equal_on_domain(partial->ma, partial->next->ma,
508
50
          partial->next->dom);
509
50
    if (
same50
)
{36
510
36
      swap_initial(sol);
511
36
      return combine_initial_into_second(sol);
512
36
    }
513
50
  }
514
763
515
690
  sol_pop_one(sol);
516
690
  sol_pop_one(sol);
517
690
518
690
  return isl_stat_ok;
519
763
}
520
521
/* Pop all solutions from the partial solution stack that were pushed onto
522
 * the stack at levels that are deeper than the current level.
523
 * If the two topmost elements on the stack have the same level
524
 * and represent the same solution, then their domains are combined.
525
 * This combined domain is the same as the current context domain
526
 * as sol_pop is called each time we move back to a higher level.
527
 * If the outer level (0) has been reached, then all partial solutions
528
 * at the current level are also popped off.
529
 */
530
static void sol_pop(struct isl_sol *sol)
531
3.77k
{
532
3.77k
  struct isl_partial_sol *partial;
533
3.77k
534
3.77k
  if (sol->error)
535
0
    return;
536
3.77k
537
3.77k
  partial = sol->partial;
538
3.77k
  if (!partial)
539
683
    return;
540
3.77k
541
3.09k
  
if (3.09k
partial->level == 0 && 3.09k
sol->level == 01.55k
)
{1.55k
542
3.11k
    for (partial = sol->partial; 
partial3.11k
;
partial = sol->partial1.55k
)
543
1.55k
      sol_pop_one(sol);
544
1.55k
    return;
545
1.55k
  }
546
3.09k
547
1.53k
  
if (1.53k
partial->level <= sol->level1.53k
)
548
4
    return;
549
1.53k
550
1.52k
  
if (1.52k
partial->next && 1.52k
partial->next->level == partial->level864
)
{772
551
772
    if (combine_initial_if_equal(sol) < 0)
552
0
      goto error;
553
772
  } else
554
756
    sol_pop_one(sol);
555
1.52k
556
1.52k
  
if (1.52k
sol->level == 01.52k
)
{646
557
747
    for (partial = sol->partial; 
partial747
;
partial = sol->partial101
)
558
101
      sol_pop_one(sol);
559
646
    return;
560
646
  }
561
1.52k
562
882
  
if (882
0882
)
563
0
error:    sol->error = 1;
564
882
}
565
566
static void sol_dec_level(struct isl_sol *sol)
567
969
{
568
969
  if (sol->error)
569
0
    return;
570
969
571
969
  sol->level--;
572
969
573
969
  sol_pop(sol);
574
969
}
575
576
static isl_stat sol_dec_level_wrap(struct isl_tab_callback *cb)
577
969
{
578
969
  struct isl_sol_callback *callback = (struct isl_sol_callback *)cb;
579
969
580
969
  sol_dec_level(callback->sol);
581
969
582
969
  return callback->sol->error ? 
isl_stat_error0
:
isl_stat_ok969
;
583
969
}
584
585
/* Move down to next level and push callback onto context tableau
586
 * to decrease the level again when it gets rolled back across
587
 * the current state.  That is, dec_level will be called with
588
 * the context tableau in the same state as it is when inc_level
589
 * is called.
590
 */
591
static void sol_inc_level(struct isl_sol *sol)
592
8.27k
{
593
8.27k
  struct isl_tab *tab;
594
8.27k
595
8.27k
  if (sol->error)
596
0
    return;
597
8.27k
598
8.27k
  sol->level++;
599
8.27k
  tab = sol->context->op->peek_tab(sol->context);
600
8.27k
  if (isl_tab_push_callback(tab, &sol->dec_level.callback) < 0)
601
0
    sol->error = 1;
602
8.27k
}
603
604
static void scale_rows(struct isl_mat *mat, isl_int m, int n_row)
605
6.23k
{
606
6.23k
  int i;
607
6.23k
608
6.23k
  if (isl_int_is_one(m))
609
6.21k
    return;
610
6.23k
611
39
  
for (i = 0; 12
i < n_row39
;
++i27
)
612
27
    isl_seq_scale(mat->row[i], mat->row[i], m, mat->n_col);
613
12
}
614
615
/* Add the solution identified by the tableau and the context tableau.
616
 *
617
 * The layout of the variables is as follows.
618
 *  tab->n_var is equal to the total number of variables in the input
619
 *      map (including divs that were copied from the context)
620
 *      + the number of extra divs constructed
621
 *      Of these, the first tab->n_param and the last tab->n_div variables
622
 *  correspond to the variables in the context, i.e.,
623
 *    tab->n_param + tab->n_div = context_tab->n_var
624
 *  tab->n_param is equal to the number of parameters and input
625
 *      dimensions in the input map
626
 *  tab->n_div is equal to the number of divs in the context
627
 *
628
 * If there is no solution, then call add_empty with a basic set
629
 * that corresponds to the context tableau.  (If add_empty is NULL,
630
 * then do nothing).
631
 *
632
 * If there is a solution, then first construct a matrix that maps
633
 * all dimensions of the context to the output variables, i.e.,
634
 * the output dimensions in the input map.
635
 * The divs in the input map (if any) that do not correspond to any
636
 * div in the context do not appear in the solution.
637
 * The algorithm will make sure that they have an integer value,
638
 * but these values themselves are of no interest.
639
 * We have to be careful not to drop or rearrange any divs in the
640
 * context because that would change the meaning of the matrix.
641
 *
642
 * To extract the value of the output variables, it should be noted
643
 * that we always use a big parameter M in the main tableau and so
644
 * the variable stored in this tableau is not an output variable x itself, but
645
 *  x' = M + x (in case of minimization)
646
 * or
647
 *  x' = M - x (in case of maximization)
648
 * If x' appears in a column, then its optimal value is zero,
649
 * which means that the optimal value of x is an unbounded number
650
 * (-M for minimization and M for maximization).
651
 * We currently assume that the output dimensions in the original map
652
 * are bounded, so this cannot occur.
653
 * Similarly, when x' appears in a row, then the coefficient of M in that
654
 * row is necessarily 1.
655
 * If the row in the tableau represents
656
 *  d x' = c + d M + e(y)
657
 * then, in case of minimization, the corresponding row in the matrix
658
 * will be
659
 *  a c + a e(y)
660
 * with a d = m, the (updated) common denominator of the matrix.
661
 * In case of maximization, the row will be
662
 *  -a c - a e(y)
663
 */
664
static void sol_add(struct isl_sol *sol, struct isl_tab *tab)
665
11.0k
{
666
11.0k
  struct isl_basic_set *bset = NULL;
667
11.0k
  struct isl_mat *mat = NULL;
668
11.0k
  unsigned off;
669
11.0k
  int row;
670
11.0k
  isl_int m;
671
11.0k
672
11.0k
  if (
sol->error || 11.0k
!tab11.0k
)
673
0
    goto error;
674
11.0k
675
11.0k
  
if (11.0k
tab->empty && 11.0k
!sol->add_empty8.12k
)
676
594
    return;
677
10.4k
  
if (10.4k
sol->context->op->is_empty(sol->context)10.4k
)
678
6.60k
    return;
679
10.4k
680
3.87k
  bset = sol_domain(sol);
681
3.87k
682
3.87k
  if (
tab->empty3.87k
)
{919
683
919
    sol_push_sol(sol, bset, NULL);
684
919
    return;
685
919
  }
686
3.87k
687
2.95k
  off = 2 + tab->M;
688
2.95k
689
2.95k
  mat = isl_mat_alloc(tab->mat->ctx, 1 + sol->n_out,
690
2.95k
              1 + tab->n_param + tab->n_div);
691
2.95k
  if (!mat)
692
0
    goto error;
693
2.95k
694
2.95k
  
isl_int_init2.95k
(m);2.95k
695
2.95k
696
2.95k
  isl_seq_clr(mat->row[0] + 1, mat->n_col - 1);
697
2.95k
  isl_int_set_si(mat->row[0][0], 1);
698
9.19k
  for (row = 0; 
row < sol->n_out9.19k
;
++row6.23k
)
{6.23k
699
6.23k
    int i = tab->n_param + row;
700
6.23k
    int r, j;
701
6.23k
702
6.23k
    isl_seq_clr(mat->row[1 + row], mat->n_col);
703
6.23k
    if (
!tab->var[i].is_row6.23k
)
{0
704
0
      if (tab->M)
705
0
        isl_die(mat->ctx, isl_error_invalid,
706
0
          "unbounded optimum", goto error2);
707
0
      continue;
708
0
    }
709
6.23k
710
6.23k
    r = tab->var[i].index;
711
6.23k
    if (tab->M &&
712
6.23k
        isl_int_ne(tab->mat->row[r][2], tab->mat->row[r][0]))
713
0
      isl_die(mat->ctx, isl_error_invalid,
714
6.23k
        "unbounded optimum", goto error2);
715
6.23k
    
isl_int_gcd6.23k
(m, mat->row[0][0], tab->mat->row[r][0]);6.23k
716
6.23k
    isl_int_divexact(m, tab->mat->row[r][0], m);
717
6.23k
    scale_rows(mat, m, 1 + row);
718
6.23k
    isl_int_divexact(m, mat->row[0][0], tab->mat->row[r][0]);
719
6.23k
    isl_int_mul(mat->row[1 + row][0], m, tab->mat->row[r][1]);
720
33.0k
    for (j = 0; 
j < tab->n_param33.0k
;
++j26.7k
)
{26.7k
721
26.7k
      int col;
722
26.7k
      if (tab->var[j].is_row)
723
14.7k
        continue;
724
12.0k
      col = tab->var[j].index;
725
12.0k
      isl_int_mul(mat->row[1 + row][1 + j], m,
726
12.0k
            tab->mat->row[r][off + col]);
727
12.0k
    }
728
7.29k
    for (j = 0; 
j < tab->n_div7.29k
;
++j1.06k
)
{1.06k
729
1.06k
      int col;
730
1.06k
      if (tab->var[tab->n_var - tab->n_div+j].is_row)
731
268
        continue;
732
796
      col = tab->var[tab->n_var - tab->n_div+j].index;
733
796
      isl_int_mul(mat->row[1 + row][1 + tab->n_param + j], m,
734
796
            tab->mat->row[r][off + col]);
735
796
    }
736
6.23k
    if (sol->max)
737
4.44k
      isl_seq_neg(mat->row[1 + row], mat->row[1 + row],
738
4.44k
            mat->n_col);
739
6.23k
  }
740
2.95k
741
2.95k
  
isl_int_clear2.95k
(m);2.95k
742
2.95k
743
2.95k
  sol_push_sol_mat(sol, bset, mat);
744
2.95k
  return;
745
0
error2:
746
0
  isl_int_clear(m);
747
0
error:
748
0
  isl_basic_set_free(bset);
749
0
  isl_mat_free(mat);
750
0
  sol->error = 1;
751
0
}
752
753
struct isl_sol_map {
754
  struct isl_sol  sol;
755
  struct isl_map  *map;
756
  struct isl_set  *empty;
757
};
758
759
static void sol_map_free(struct isl_sol *sol)
760
1.27k
{
761
1.27k
  struct isl_sol_map *sol_map = (struct isl_sol_map *) sol;
762
1.27k
  isl_map_free(sol_map->map);
763
1.27k
  isl_set_free(sol_map->empty);
764
1.27k
}
765
766
/* This function is called for parts of the context where there is
767
 * no solution, with "bset" corresponding to the context tableau.
768
 * Simply add the basic set to the set "empty".
769
 */
770
static void sol_map_add_empty(struct isl_sol_map *sol,
771
  struct isl_basic_set *bset)
772
766
{
773
766
  if (
!bset || 766
!sol->empty766
)
774
0
    goto error;
775
766
776
766
  sol->empty = isl_set_grow(sol->empty, 1);
777
766
  bset = isl_basic_set_simplify(bset);
778
766
  bset = isl_basic_set_finalize(bset);
779
766
  sol->empty = isl_set_add_basic_set(sol->empty, isl_basic_set_copy(bset));
780
766
  if (!sol->empty)
781
0
    goto error;
782
766
  isl_basic_set_free(bset);
783
766
  return;
784
0
error:
785
0
  isl_basic_set_free(bset);
786
0
  sol->sol.error = 1;
787
0
}
788
789
static void sol_map_add_empty_wrap(struct isl_sol *sol,
790
  struct isl_basic_set *bset)
791
766
{
792
766
  sol_map_add_empty((struct isl_sol_map *)sol, bset);
793
766
}
794
795
/* Given a basic set "dom" that represents the context and a tuple of
796
 * affine expressions "ma" defined over this domain, construct a basic map
797
 * that expresses this function on the domain.
798
 */
799
static void sol_map_add(struct isl_sol_map *sol,
800
  __isl_take isl_basic_set *dom, __isl_take isl_multi_aff *ma)
801
1.21k
{
802
1.21k
  isl_basic_map *bmap;
803
1.21k
804
1.21k
  if (
sol->sol.error || 1.21k
!dom1.21k
||
!ma1.21k
)
805
0
    goto error;
806
1.21k
807
1.21k
  bmap = isl_basic_map_from_multi_aff2(ma, sol->sol.rational);
808
1.21k
  bmap = isl_basic_map_intersect_domain(bmap, dom);
809
1.21k
  sol->map = isl_map_grow(sol->map, 1);
810
1.21k
  sol->map = isl_map_add_basic_map(sol->map, bmap);
811
1.21k
  if (!sol->map)
812
0
    sol->sol.error = 1;
813
1.21k
  return;
814
0
error:
815
0
  isl_basic_set_free(dom);
816
0
  isl_multi_aff_free(ma);
817
0
  sol->sol.error = 1;
818
0
}
819
820
static void sol_map_add_wrap(struct isl_sol *sol,
821
  __isl_take isl_basic_set *dom, __isl_take isl_multi_aff *ma)
822
1.21k
{
823
1.21k
  sol_map_add((struct isl_sol_map *)sol, dom, ma);
824
1.21k
}
825
826
827
/* Store the "parametric constant" of row "row" of tableau "tab" in "line",
828
 * i.e., the constant term and the coefficients of all variables that
829
 * appear in the context tableau.
830
 * Note that the coefficient of the big parameter M is NOT copied.
831
 * The context tableau may not have a big parameter and even when it
832
 * does, it is a different big parameter.
833
 */
834
static void get_row_parameter_line(struct isl_tab *tab, int row, isl_int *line)
835
12.1k
{
836
12.1k
  int i;
837
12.1k
  unsigned off = 2 + tab->M;
838
12.1k
839
12.1k
  isl_int_set(line[0], tab->mat->row[row][1]);
840
69.5k
  for (i = 0; 
i < tab->n_param69.5k
;
++i57.4k
)
{57.4k
841
57.4k
    if (tab->var[i].is_row)
842
27.3k
      isl_int_set_si(line[1 + i], 0);
843
30.0k
    else {
844
30.0k
      int col = tab->var[i].index;
845
30.0k
      isl_int_set(line[1 + i], tab->mat->row[row][off + col]);
846
30.0k
    }
847
57.4k
  }
848
20.3k
  for (i = 0; 
i < tab->n_div20.3k
;
++i8.24k
)
{8.24k
849
8.24k
    if (tab->var[tab->n_var - tab->n_div + i].is_row)
850
3.26k
      isl_int_set_si(line[1 + tab->n_param + i], 0);
851
4.98k
    else {
852
4.98k
      int col = tab->var[tab->n_var - tab->n_div + i].index;
853
4.98k
      isl_int_set(line[1 + tab->n_param + i],
854
4.98k
            tab->mat->row[row][off + col]);
855
4.98k
    }
856
8.24k
  }
857
12.1k
}
858
859
/* Check if rows "row1" and "row2" have identical "parametric constants",
860
 * as explained above.
861
 * In this case, we also insist that the coefficients of the big parameter
862
 * be the same as the values of the constants will only be the same
863
 * if these coefficients are also the same.
864
 */
865
static int identical_parameter_line(struct isl_tab *tab, int row1, int row2)
866
38.8k
{
867
38.8k
  int i;
868
38.8k
  unsigned off = 2 + tab->M;
869
38.8k
870
38.8k
  if (isl_int_ne(tab->mat->row[row1][1], tab->mat->row[row2][1]))
871
35.6k
    return 0;
872
38.8k
873
3.20k
  
if (3.20k
tab->M && 3.20k
isl_int_ne3.20k
(tab->mat->row[row1][2],
874
3.20k
         tab->mat->row[row2][2]))
875
1.07k
    return 0;
876
3.20k
877
4.44k
  
for (i = 0; 2.12k
i < tab->n_param + tab->n_div4.44k
;
++i2.31k
)
{4.39k
878
4.24k
    int pos = i < tab->n_param ? i :
879
153
      tab->n_var - tab->n_div + i - tab->n_param;
880
4.39k
    int col;
881
4.39k
882
4.39k
    if (tab->var[pos].is_row)
883
912
      continue;
884
3.48k
    col = tab->var[pos].index;
885
3.48k
    if (isl_int_ne(tab->mat->row[row1][off + col],
886
3.48k
             tab->mat->row[row2][off + col]))
887
2.08k
      return 0;
888
3.48k
  }
889
47
  return 1;
890
2.12k
}
891
892
/* Return an inequality that expresses that the "parametric constant"
893
 * should be non-negative.
894
 * This function is only called when the coefficient of the big parameter
895
 * is equal to zero.
896
 */
897
static struct isl_vec *get_row_parameter_ineq(struct isl_tab *tab, int row)
898
8.17k
{
899
8.17k
  struct isl_vec *ineq;
900
8.17k
901
8.17k
  ineq = isl_vec_alloc(tab->mat->ctx, 1 + tab->n_param + tab->n_div);
902
8.17k
  if (!ineq)
903
0
    return NULL;
904
8.17k
905
8.17k
  get_row_parameter_line(tab, row, ineq->el);
906
8.17k
  if (ineq)
907
8.17k
    ineq = isl_vec_normalize(ineq);
908
8.17k
909
8.17k
  return ineq;
910
8.17k
}
911
912
/* Normalize a div expression of the form
913
 *
914
 *  [(g*f(x) + c)/(g * m)]
915
 *
916
 * with c the constant term and f(x) the remaining coefficients, to
917
 *
918
 *  [(f(x) + [c/g])/m]
919
 */
920
static void normalize_div(__isl_keep isl_vec *div)
921
401
{
922
401
  isl_ctx *ctx = isl_vec_get_ctx(div);
923
401
  int len = div->size - 2;
924
401
925
401
  isl_seq_gcd(div->el + 2, len, &ctx->normalize_gcd);
926
401
  isl_int_gcd(ctx->normalize_gcd, ctx->normalize_gcd, div->el[0]);
927
401
928
401
  if (isl_int_is_one(ctx->normalize_gcd))
929
383
    return;
930
401
931
18
  
isl_int_divexact18
(div->el[0], div->el[0], ctx->normalize_gcd);18
932
18
  isl_int_fdiv_q(div->el[1], div->el[1], ctx->normalize_gcd);
933
18
  isl_seq_scale_down(div->el + 2, div->el + 2, ctx->normalize_gcd, len);
934
18
}
935
936
/* Return an integer division for use in a parametric cut based
937
 * on the given row.
938
 * In particular, let the parametric constant of the row be
939
 *
940
 *    \sum_i a_i y_i
941
 *
942
 * where y_0 = 1, but none of the y_i corresponds to the big parameter M.
943
 * The div returned is equal to
944
 *
945
 *    floor(\sum_i {-a_i} y_i) = floor((\sum_i (-a_i mod d) y_i)/d)
946
 */
947
static struct isl_vec *get_row_parameter_div(struct isl_tab *tab, int row)
948
365
{
949
365
  struct isl_vec *div;
950
365
951
365
  div = isl_vec_alloc(tab->mat->ctx, 1 + 1 + tab->n_param + tab->n_div);
952
365
  if (!div)
953
0
    return NULL;
954
365
955
365
  
isl_int_set365
(div->el[0], tab->mat->row[row][0]);365
956
365
  get_row_parameter_line(tab, row, div->el + 1);
957
365
  isl_seq_neg(div->el + 1, div->el + 1, div->size - 1);
958
365
  normalize_div(div);
959
365
  isl_seq_fdiv_r(div->el + 1, div->el + 1, div->el[0], div->size - 1);
960
365
961
365
  return div;
962
365
}
963
964
/* Return an integer division for use in transferring an integrality constraint
965
 * to the context.
966
 * In particular, let the parametric constant of the row be
967
 *
968
 *    \sum_i a_i y_i
969
 *
970
 * where y_0 = 1, but none of the y_i corresponds to the big parameter M.
971
 * The the returned div is equal to
972
 *
973
 *    floor(\sum_i {a_i} y_i) = floor((\sum_i (a_i mod d) y_i)/d)
974
 */
975
static struct isl_vec *get_row_split_div(struct isl_tab *tab, int row)
976
36
{
977
36
  struct isl_vec *div;
978
36
979
36
  div = isl_vec_alloc(tab->mat->ctx, 1 + 1 + tab->n_param + tab->n_div);
980
36
  if (!div)
981
0
    return NULL;
982
36
983
36
  
isl_int_set36
(div->el[0], tab->mat->row[row][0]);36
984
36
  get_row_parameter_line(tab, row, div->el + 1);
985
36
  normalize_div(div);
986
36
  isl_seq_fdiv_r(div->el + 1, div->el + 1, div->el[0], div->size - 1);
987
36
988
36
  return div;
989
36
}
990
991
/* Construct and return an inequality that expresses an upper bound
992
 * on the given div.
993
 * In particular, if the div is given by
994
 *
995
 *  d = floor(e/m)
996
 *
997
 * then the inequality expresses
998
 *
999
 *  m d <= e
1000
 */
1001
static struct isl_vec *ineq_for_div(struct isl_basic_set *bset, unsigned div)
1002
36
{
1003
36
  unsigned total;
1004
36
  unsigned div_pos;
1005
36
  struct isl_vec *ineq;
1006
36
1007
36
  if (!bset)
1008
0
    return NULL;
1009
36
1010
36
  total = isl_basic_set_total_dim(bset);
1011
36
  div_pos = 1 + total - bset->n_div + div;
1012
36
1013
36
  ineq = isl_vec_alloc(bset->ctx, 1 + total);
1014
36
  if (!ineq)
1015
0
    return NULL;
1016
36
1017
36
  isl_seq_cpy(ineq->el, bset->div[div] + 1, 1 + total);
1018
36
  isl_int_neg(ineq->el[div_pos], bset->div[div][0]);
1019
36
  return ineq;
1020
36
}
1021
1022
/* Given a row in the tableau and a div that was created
1023
 * using get_row_split_div and that has been constrained to equality, i.e.,
1024
 *
1025
 *    d = floor(\sum_i {a_i} y_i) = \sum_i {a_i} y_i
1026
 *
1027
 * replace the expression "\sum_i {a_i} y_i" in the row by d,
1028
 * i.e., we subtract "\sum_i {a_i} y_i" and add 1 d.
1029
 * The coefficients of the non-parameters in the tableau have been
1030
 * verified to be integral.  We can therefore simply replace coefficient b
1031
 * by floor(b).  For the coefficients of the parameters we have
1032
 * floor(a_i) = a_i - {a_i}, while for the other coefficients, we have
1033
 * floor(b) = b.
1034
 */
1035
static struct isl_tab *set_row_cst_to_div(struct isl_tab *tab, int row, int div)
1036
36
{
1037
36
  isl_seq_fdiv_q(tab->mat->row[row] + 1, tab->mat->row[row] + 1,
1038
36
      tab->mat->row[row][0], 1 + tab->M + tab->n_col);
1039
36
1040
36
  isl_int_set_si(tab->mat->row[row][0], 1);
1041
36
1042
36
  if (
tab->var[tab->n_var - tab->n_div + div].is_row36
)
{0
1043
0
    int drow = tab->var[tab->n_var - tab->n_div + div].index;
1044
0
1045
0
    isl_assert(tab->mat->ctx,
1046
0
      isl_int_is_one(tab->mat->row[drow][0]), goto error);
1047
0
    isl_seq_combine(tab->mat->row[row] + 1,
1048
0
      tab->mat->ctx->one, tab->mat->row[row] + 1,
1049
0
      tab->mat->ctx->one, tab->mat->row[drow] + 1,
1050
0
      1 + tab->M + tab->n_col);
1051
36
  } else {
1052
36
    int dcol = tab->var[tab->n_var - tab->n_div + div].index;
1053
36
1054
36
    isl_int_add_ui(tab->mat->row[row][2 + tab->M + dcol],
1055
36
        tab->mat->row[row][2 + tab->M + dcol], 1);
1056
36
  }
1057
36
1058
36
  return tab;
1059
0
error:
1060
0
  isl_tab_free(tab);
1061
0
  return NULL;
1062
36
}
1063
1064
/* Check if the (parametric) constant of the given row is obviously
1065
 * negative, meaning that we don't need to consult the context tableau.
1066
 * If there is a big parameter and its coefficient is non-zero,
1067
 * then this coefficient determines the outcome.
1068
 * Otherwise, we check whether the constant is negative and
1069
 * all non-zero coefficients of parameters are negative and
1070
 * belong to non-negative parameters.
1071
 */
1072
static int is_obviously_neg(struct isl_tab *tab, int row)
1073
147k
{
1074
147k
  int i;
1075
147k
  int col;
1076
147k
  unsigned off = 2 + tab->M;
1077
147k
1078
147k
  if (
tab->M147k
)
{38.1k
1079
38.1k
    if (isl_int_is_pos(tab->mat->row[row][2]))
1080
19.2k
      return 0;
1081
18.9k
    
if (18.9k
isl_int_is_neg18.9k
(tab->mat->row[row][2]))
1082
0
      return 1;
1083
18.9k
  }
1084
147k
1085
128k
  
if (128k
isl_int_is_nonneg128k
(tab->mat->row[row][1]))
1086
119k
    return 0;
1087
15.5k
  
for (i = 0; 8.70k
i < tab->n_param15.5k
;
++i6.80k
)
{11.6k
1088
11.6k
    /* Eliminated parameter */
1089
11.6k
    if (tab->var[i].is_row)
1090
3.54k
      continue;
1091
8.05k
    col = tab->var[i].index;
1092
8.05k
    if (isl_int_is_zero(tab->mat->row[row][off + col]))
1093
3.14k
      continue;
1094
4.90k
    
if (4.90k
!tab->var[i].is_nonneg4.90k
)
1095
3.87k
      return 0;
1096
1.03k
    
if (1.03k
isl_int_is_pos1.03k
(tab->mat->row[row][off + col]))
1097
931
      return 0;
1098
1.03k
  }
1099
4.32k
  
for (i = 0; 3.89k
i < tab->n_div4.32k
;
++i424
)
{493
1100
493
    if (tab->var[tab->n_var - tab->n_div + i].is_row)
1101
403
      continue;
1102
90
    col = tab->var[tab->n_var - tab->n_div + i].index;
1103
90
    if (isl_int_is_zero(tab->mat->row[row][off + col]))
1104
21
      continue;
1105
69
    
if (69
!tab->var[tab->n_var - tab->n_div + i].is_nonneg69
)
1106
55
      return 0;
1107
14
    
if (14
isl_int_is_pos14
(tab->mat->row[row][off + col]))
1108
14
      return 0;
1109
14
  }
1110
3.82k
  return 1;
1111
3.89k
}
1112
1113
/* Check if the (parametric) constant of the given row is obviously
1114
 * non-negative, meaning that we don't need to consult the context tableau.
1115
 * If there is a big parameter and its coefficient is non-zero,
1116
 * then this coefficient determines the outcome.
1117
 * Otherwise, we check whether the constant is non-negative and
1118
 * all non-zero coefficients of parameters are positive and
1119
 * belong to non-negative parameters.
1120
 */
1121
static int is_obviously_nonneg(struct isl_tab *tab, int row)
1122
14.7k
{
1123
14.7k
  int i;
1124
14.7k
  int col;
1125
14.7k
  unsigned off = 2 + tab->M;
1126
14.7k
1127
14.7k
  if (
tab->M14.7k
)
{14.7k
1128
14.7k
    if (isl_int_is_pos(tab->mat->row[row][2]))
1129
5.79k
      return 1;
1130
8.95k
    
if (8.95k
isl_int_is_neg8.95k
(tab->mat->row[row][2]))
1131
0
      return 0;
1132
8.95k
  }
1133
14.7k
1134
8.95k
  
if (8.95k
isl_int_is_neg8.95k
(tab->mat->row[row][1]))
1135
2.95k
    return 0;
1136
21.8k
  
for (i = 0; 5.99k
i < tab->n_param21.8k
;
++i15.8k
)
{18.2k
1137
18.2k
    /* Eliminated parameter */
1138
18.2k
    if (tab->var[i].is_row)
1139
7.33k
      continue;
1140
10.9k
    col = tab->var[i].index;
1141
10.9k
    if (isl_int_is_zero(tab->mat->row[row][off + col]))
1142
7.41k
      continue;
1143
3.53k
    
if (3.53k
!tab->var[i].is_nonneg3.53k
)
1144
798
      return 0;
1145
2.73k
    
if (2.73k
isl_int_is_neg2.73k
(tab->mat->row[row][off + col]))
1146
1.65k
      return 0;
1147
2.73k
  }
1148
8.18k
  
for (i = 0; 3.54k
i < tab->n_div8.18k
;
++i4.64k
)
{4.71k
1149
4.71k
    if (tab->var[tab->n_var - tab->n_div + i].is_row)
1150
3.89k
      continue;
1151
823
    col = tab->var[tab->n_var - tab->n_div + i].index;
1152
823
    if (isl_int_is_zero(tab->mat->row[row][off + col]))
1153
749
      continue;
1154
74
    
if (74
!tab->var[tab->n_var - tab->n_div + i].is_nonneg74
)
1155
68
      return 0;
1156
6
    
if (6
isl_int_is_neg6
(tab->mat->row[row][off + col]))
1157
3
      return 0;
1158
6
  }
1159
3.47k
  return 1;
1160
3.54k
}
1161
1162
/* Given a row r and two columns, return the column that would
1163
 * lead to the lexicographically smallest increment in the sample
1164
 * solution when leaving the basis in favor of the row.
1165
 * Pivoting with column c will increment the sample value by a non-negative
1166
 * constant times a_{V,c}/a_{r,c}, with a_{V,c} the elements of column c
1167
 * corresponding to the non-parametric variables.
1168
 * If variable v appears in a column c_v, the a_{v,c} = 1 iff c = c_v,
1169
 * with all other entries in this virtual row equal to zero.
1170
 * If variable v appears in a row, then a_{v,c} is the element in column c
1171
 * of that row.
1172
 *
1173
 * Let v be the first variable with a_{v,c1}/a_{r,c1} != a_{v,c2}/a_{r,c2}.
1174
 * Then if a_{v,c1}/a_{r,c1} < a_{v,c2}/a_{r,c2}, i.e.,
1175
 * a_{v,c2} a_{r,c1} - a_{v,c1} a_{r,c2} > 0, c1 results in the minimal
1176
 * increment.  Otherwise, it's c2.
1177
 */
1178
static int lexmin_col_pair(struct isl_tab *tab,
1179
  int row, int col1, int col2, isl_int tmp)
1180
6.45k
{
1181
6.45k
  int i;
1182
6.45k
  isl_int *tr;
1183
6.45k
1184
6.45k
  tr = tab->mat->row[row] + 2 + tab->M;
1185
6.45k
1186
32.1k
  for (i = tab->n_param; 
i < tab->n_var - tab->n_div32.1k
;
++i25.6k
)
{32.1k
1187
32.1k
    int s1, s2;
1188
32.1k
    isl_int *r;
1189
32.1k
1190
32.1k
    if (
!tab->var[i].is_row32.1k
)
{11.9k
1191
11.9k
      if (tab->var[i].index == col1)
1192
827
        return col2;
1193
11.1k
      
if (11.1k
tab->var[i].index == col211.1k
)
1194
898
        return col1;
1195
10.2k
      continue;
1196
11.1k
    }
1197
32.1k
1198
20.1k
    
if (20.1k
tab->var[i].index == row20.1k
)
1199
173
      continue;
1200
20.1k
1201
19.9k
    r = tab->mat->row[tab->var[i].index] + 2 + tab->M;
1202
19.9k
    s1 = isl_int_sgn(r[col1]);
1203
19.9k
    s2 = isl_int_sgn(r[col2]);
1204
19.9k
    if (
s1 == 0 && 19.9k
s2 == 016.4k
)
1205
14.1k
      continue;
1206
5.86k
    
if (5.86k
s1 < s25.86k
)
1207
2.53k
      return col1;
1208
3.33k
    
if (3.33k
s2 < s13.33k
)
1209
776
      return col2;
1210
3.33k
1211
2.55k
    
isl_int_mul2.55k
(tmp, r[col2], tr[col1]);2.55k
1212
2.55k
    isl_int_submul(tmp, r[col1], tr[col2]);
1213
2.55k
    if (isl_int_is_pos(tmp))
1214
901
      return col1;
1215
1.65k
    
if (1.65k
isl_int_is_neg1.65k
(tmp))
1216
518
      return col2;
1217
1.65k
  }
1218
0
  return -1;
1219
6.45k
}
1220
1221
/* Given a row in the tableau, find and return the column that would
1222
 * result in the lexicographically smallest, but positive, increment
1223
 * in the sample point.
1224
 * If there is no such column, then return tab->n_col.
1225
 * If anything goes wrong, return -1.
1226
 */
1227
static int lexmin_pivot_col(struct isl_tab *tab, int row)
1228
7.47k
{
1229
7.47k
  int j;
1230
7.47k
  int col = tab->n_col;
1231
7.47k
  isl_int *tr;
1232
7.47k
  isl_int tmp;
1233
7.47k
1234
7.47k
  tr = tab->mat->row[row] + 2 + tab->M;
1235
7.47k
1236
7.47k
  isl_int_init(tmp);
1237
7.47k
1238
72.1k
  for (j = tab->n_dead; 
j < tab->n_col72.1k
;
++j64.6k
)
{64.6k
1239
64.6k
    if (tab->col_var[j] >= 0 &&
1240
46.5k
        (tab->col_var[j] < tab->n_param  ||
1241
39.0k
        tab->col_var[j] >= tab->n_var - tab->n_div))
1242
9.42k
      continue;
1243
64.6k
1244
55.2k
    
if (55.2k
!55.2k
isl_int_is_pos55.2k
(tr[j]))
1245
42.7k
      continue;
1246
55.2k
1247
12.4k
    
if (12.4k
col == tab->n_col12.4k
)
1248
6.04k
      col = j;
1249
12.4k
    else
1250
6.45k
      col = lexmin_col_pair(tab, row, col, j, tmp);
1251
12.4k
    isl_assert(tab->mat->ctx, col >= 0, goto error);
1252
12.4k
  }
1253
7.47k
1254
7.47k
  
isl_int_clear7.47k
(tmp);7.47k
1255
7.47k
  return col;
1256
0
error:
1257
0
  isl_int_clear(tmp);
1258
0
  return -1;
1259
7.47k
}
1260
1261
/* Return the first known violated constraint, i.e., a non-negative
1262
 * constraint that currently has an either obviously negative value
1263
 * or a previously determined to be negative value.
1264
 *
1265
 * If any constraint has a negative coefficient for the big parameter,
1266
 * if any, then we return one of these first.
1267
 */
1268
static int first_neg(struct isl_tab *tab)
1269
21.6k
{
1270
21.6k
  int row;
1271
21.6k
1272
21.6k
  if (tab->M)
1273
116k
    
for (row = tab->n_redundant; 13.5k
row < tab->n_row116k
;
++row103k
)
{104k
1274
104k
      if (!isl_tab_var_from_row(tab, row)->is_nonneg)
1275
22.5k
        continue;
1276
82.4k
      
if (82.4k
!82.4k
isl_int_is_neg82.4k
(tab->mat->row[row][2]))
1277
80.7k
        continue;
1278
1.71k
      
if (1.71k
tab->row_sign1.71k
)
1279
1.71k
        tab->row_sign[row] = isl_tab_row_neg;
1280
1.71k
      return row;
1281
82.4k
    }
1282
216k
  
for (row = tab->n_redundant; 19.9k
row < tab->n_row216k
;
++row196k
)
{202k
1283
202k
    if (!isl_tab_var_from_row(tab, row)->is_nonneg)
1284
19.5k
      continue;
1285
183k
    
if (183k
tab->row_sign183k
)
{73.5k
1286
73.5k
      if (tab->row_sign[row] == 0 &&
1287
38.1k
          is_obviously_neg(tab, row))
1288
203
        tab->row_sign[row] = isl_tab_row_neg;
1289
73.5k
      if (tab->row_sign[row] != isl_tab_row_neg)
1290
71.4k
        continue;
1291
109k
    } else 
if (109k
!is_obviously_neg(tab, row)109k
)
1292
105k
      continue;
1293
5.76k
    return row;
1294
183k
  }
1295
14.1k
  return -1;
1296
19.9k
}
1297
1298
/* Check whether the invariant that all columns are lexico-positive
1299
 * is satisfied.  This function is not called from the current code
1300
 * but is useful during debugging.
1301
 */
1302
static void check_lexpos(struct isl_tab *tab) __attribute__ ((unused));
1303
static void check_lexpos(struct isl_tab *tab)
1304
0
{
1305
0
  unsigned off = 2 + tab->M;
1306
0
  int col;
1307
0
  int var;
1308
0
  int row;
1309
0
1310
0
  for (col = tab->n_dead; col < tab->n_col; ++col) {
1311
0
    if (tab->col_var[col] >= 0 &&
1312
0
        (tab->col_var[col] < tab->n_param ||
1313
0
         tab->col_var[col] >= tab->n_var - tab->n_div))
1314
0
      continue;
1315
0
    for (var = tab->n_param; var < tab->n_var - tab->n_div; ++var) {
1316
0
      if (!tab->var[var].is_row) {
1317
0
        if (tab->var[var].index == col)
1318
0
          break;
1319
0
        else
1320
0
          continue;
1321
0
      }
1322
0
      row = tab->var[var].index;
1323
0
      if (isl_int_is_zero(tab->mat->row[row][off + col]))
1324
0
        continue;
1325
0
      if (isl_int_is_pos(tab->mat->row[row][off + col]))
1326
0
        break;
1327
0
      fprintf(stderr, "lexneg column %d (row %d)\n",
1328
0
        col, row);
1329
0
    }
1330
0
    if (var >= tab->n_var - tab->n_div)
1331
0
      fprintf(stderr, "zero column %d\n", col);
1332
0
  }
1333
0
}
1334
1335
/* Report to the caller that the given constraint is part of an encountered
1336
 * conflict.
1337
 */
1338
static int report_conflicting_constraint(struct isl_tab *tab, int con)
1339
1.17k
{
1340
1.17k
  return tab->conflict(con, tab->conflict_user);
1341
1.17k
}
1342
1343
/* Given a conflicting row in the tableau, report all constraints
1344
 * involved in the row to the caller.  That is, the row itself
1345
 * (if it represents a constraint) and all constraint columns with
1346
 * non-zero (and therefore negative) coefficients.
1347
 */
1348
static int report_conflict(struct isl_tab *tab, int row)
1349
1.43k
{
1350
1.43k
  int j;
1351
1.43k
  isl_int *tr;
1352
1.43k
1353
1.43k
  if (!tab->conflict)
1354
995
    return 0;
1355
1.43k
1356
436
  
if (436
tab->row_var[row] < 0 &&436
1357
436
      report_conflicting_constraint(tab, ~tab->row_var[row]) < 0)
1358
0
    return -1;
1359
436
1360
436
  tr = tab->mat->row[row] + 2 + tab->M;
1361
436
1362
5.17k
  for (j = tab->n_dead; 
j < tab->n_col5.17k
;
++j4.73k
)
{4.73k
1363
4.73k
    if (tab->col_var[j] >= 0 &&
1364
3.80k
        (tab->col_var[j] < tab->n_param  ||
1365
3.80k
        tab->col_var[j] >= tab->n_var - tab->n_div))
1366
0
      continue;
1367
4.73k
1368
4.73k
    
if (4.73k
!4.73k
isl_int_is_neg4.73k
(tr[j]))
1369
3.93k
      continue;
1370
4.73k
1371
798
    
if (798
tab->col_var[j] < 0 &&798
1372
737
        report_conflicting_constraint(tab, ~tab->col_var[j]) < 0)
1373
0
      return -1;
1374
798
  }
1375
436
1376
436
  return 0;
1377
436
}
1378
1379
/* Resolve all known or obviously violated constraints through pivoting.
1380
 * In particular, as long as we can find any violated constraint, we
1381
 * look for a pivoting column that would result in the lexicographically
1382
 * smallest increment in the sample point.  If there is no such column
1383
 * then the tableau is infeasible.
1384
 */
1385
static int restore_lexmin(struct isl_tab *tab) WARN_UNUSED;
1386
static int restore_lexmin(struct isl_tab *tab)
1387
15.6k
{
1388
15.6k
  int row, col;
1389
15.6k
1390
15.6k
  if (!tab)
1391
0
    return -1;
1392
15.6k
  
if (15.6k
tab->empty15.6k
)
1393
32
    return 0;
1394
21.6k
  
while (15.6k
(row = first_neg(tab)) != -121.6k
)
{7.47k
1395
7.47k
    col = lexmin_pivot_col(tab, row);
1396
7.47k
    if (
col >= tab->n_col7.47k
)
{1.43k
1397
1.43k
      if (report_conflict(tab, row) < 0)
1398
0
        return -1;
1399
1.43k
      
if (1.43k
isl_tab_mark_empty(tab) < 01.43k
)
1400
0
        return -1;
1401
1.43k
      return 0;
1402
1.43k
    }
1403
6.04k
    
if (6.04k
col < 06.04k
)
1404
0
      return -1;
1405
6.04k
    
if (6.04k
isl_tab_pivot(tab, row, col) < 06.04k
)
1406
0
      return -1;
1407
6.04k
  }
1408
14.1k
  return 0;
1409
15.6k
}
1410
1411
/* Given a row that represents an equality, look for an appropriate
1412
 * pivoting column.
1413
 * In particular, if there are any non-zero coefficients among
1414
 * the non-parameter variables, then we take the last of these
1415
 * variables.  Eliminating this variable in terms of the other
1416
 * variables and/or parameters does not influence the property
1417
 * that all column in the initial tableau are lexicographically
1418
 * positive.  The row corresponding to the eliminated variable
1419
 * will only have non-zero entries below the diagonal of the
1420
 * initial tableau.  That is, we transform
1421
 *
1422
 *    I       I
1423
 *      1   into    a
1424
 *        I         I
1425
 *
1426
 * If there is no such non-parameter variable, then we are dealing with
1427
 * pure parameter equality and we pick any parameter with coefficient 1 or -1
1428
 * for elimination.  This will ensure that the eliminated parameter
1429
 * always has an integer value whenever all the other parameters are integral.
1430
 * If there is no such parameter then we return -1.
1431
 */
1432
static int last_var_col_or_int_par_col(struct isl_tab *tab, int row)
1433
9.47k
{
1434
9.47k
  unsigned off = 2 + tab->M;
1435
9.47k
  int i;
1436
9.47k
1437
38.5k
  for (i = tab->n_var - tab->n_div - 1; 
i >= 0 && 38.5k
i >= tab->n_param37.9k
;
--i29.0k
)
{34.4k
1438
34.4k
    int col;
1439
34.4k
    if (tab->var[i].is_row)
1440
17.8k
      continue;
1441
16.6k
    col = tab->var[i].index;
1442
16.6k
    if (col <= tab->n_dead)
1443
1.30k
      continue;
1444
15.2k
    
if (15.2k
!15.2k
isl_int_is_zero15.2k
(tab->mat->row[row][off + col]))
1445
5.35k
      return col;
1446
15.2k
  }
1447
13.1k
  
for (i = tab->n_dead; 4.12k
i < tab->n_col13.1k
;
++i9.04k
)
{13.0k
1448
13.0k
    if (isl_int_is_one(tab->mat->row[row][off + i]))
1449
1.85k
      return i;
1450
11.1k
    
if (11.1k
isl_int_is_negone11.1k
(tab->mat->row[row][off + i]))
1451
2.10k
      return i;
1452
11.1k
  }
1453
156
  return -1;
1454
4.12k
}
1455
1456
/* Add an equality that is known to be valid to the tableau.
1457
 * We first check if we can eliminate a variable or a parameter.
1458
 * If not, we add the equality as two inequalities.
1459
 * In this case, the equality was a pure parameter equality and there
1460
 * is no need to resolve any constraint violations.
1461
 *
1462
 * This function assumes that at least two more rows and at least
1463
 * two more elements in the constraint array are available in the tableau.
1464
 */
1465
static struct isl_tab *add_lexmin_valid_eq(struct isl_tab *tab, isl_int *eq)
1466
9.47k
{
1467
9.47k
  int i;
1468
9.47k
  int r;
1469
9.47k
1470
9.47k
  if (!tab)
1471
0
    return NULL;
1472
9.47k
  r = isl_tab_add_row(tab, eq);
1473
9.47k
  if (r < 0)
1474
0
    goto error;
1475
9.47k
1476
9.47k
  r = tab->con[r].index;
1477
9.47k
  i = last_var_col_or_int_par_col(tab, r);
1478
9.47k
  if (
i < 09.47k
)
{156
1479
156
    tab->con[r].is_nonneg = 1;
1480
156
    if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r]) < 0)
1481
0
      goto error;
1482
156
    isl_seq_neg(eq, eq, 1 + tab->n_var);
1483
156
    r = isl_tab_add_row(tab, eq);
1484
156
    if (r < 0)
1485
0
      goto error;
1486
156
    tab->con[r].is_nonneg = 1;
1487
156
    if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r]) < 0)
1488
0
      goto error;
1489
9.32k
  } else {
1490
9.32k
    if (isl_tab_pivot(tab, r, i) < 0)
1491
0
      goto error;
1492
9.32k
    
if (9.32k
isl_tab_kill_col(tab, i) < 09.32k
)
1493
0
      goto error;
1494
9.32k
    tab->n_eq++;
1495
9.32k
  }
1496
9.47k
1497
9.47k
  return tab;
1498
0
error:
1499
0
  isl_tab_free(tab);
1500
0
  return NULL;
1501
9.47k
}
1502
1503
/* Check if the given row is a pure constant.
1504
 */
1505
static int is_constant(struct isl_tab *tab, int row)
1506
353
{
1507
353
  unsigned off = 2 + tab->M;
1508
353
1509
353
  return isl_seq_first_non_zero(tab->mat->row[row] + off + tab->n_dead,
1510
353
          tab->n_col - tab->n_dead) == -1;
1511
353
}
1512
1513
/* Add an equality that may or may not be valid to the tableau.
1514
 * If the resulting row is a pure constant, then it must be zero.
1515
 * Otherwise, the resulting tableau is empty.
1516
 *
1517
 * If the row is not a pure constant, then we add two inequalities,
1518
 * each time checking that they can be satisfied.
1519
 * In the end we try to use one of the two constraints to eliminate
1520
 * a column.
1521
 *
1522
 * This function assumes that at least two more rows and at least
1523
 * two more elements in the constraint array are available in the tableau.
1524
 */
1525
static int add_lexmin_eq(struct isl_tab *tab, isl_int *eq) WARN_UNUSED;
1526
static int add_lexmin_eq(struct isl_tab *tab, isl_int *eq)
1527
353
{
1528
353
  int r1, r2;
1529
353
  int row;
1530
353
  struct isl_tab_undo *snap;
1531
353
1532
353
  if (!tab)
1533
0
    return -1;
1534
353
  snap = isl_tab_snap(tab);
1535
353
  r1 = isl_tab_add_row(tab, eq);
1536
353
  if (r1 < 0)
1537
0
    return -1;
1538
353
  tab->con[r1].is_nonneg = 1;
1539
353
  if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r1]) < 0)
1540
0
    return -1;
1541
353
1542
353
  row = tab->con[r1].index;
1543
353
  if (
is_constant(tab, row)353
)
{51
1544
51
    if (
!51
isl_int_is_zero51
(tab->mat->row[row][1]) ||
1545
51
        
(tab->M && 51
!0
isl_int_is_zero0
(tab->mat->row[row][2])))
{0
1546
0
      if (isl_tab_mark_empty(tab) < 0)
1547
0
        return -1;
1548
0
      return 0;
1549
0
    }
1550
51
    
if (51
isl_tab_rollback(tab, snap) < 051
)
1551
0
      return -1;
1552
51
    return 0;
1553
51
  }
1554
353
1555
302
  
if (302
restore_lexmin(tab) < 0302
)
1556
0
    return -1;
1557
302
  
if (302
tab->empty302
)
1558
20
    return 0;
1559
302
1560
282
  isl_seq_neg(eq, eq, 1 + tab->n_var);
1561
282
1562
282
  r2 = isl_tab_add_row(tab, eq);
1563
282
  if (r2 < 0)
1564
0
    return -1;
1565
282
  tab->con[r2].is_nonneg = 1;
1566
282
  if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r2]) < 0)
1567
0
    return -1;
1568
282
1569
282
  
if (282
restore_lexmin(tab) < 0282
)
1570
0
    return -1;
1571
282
  
if (282
tab->empty282
)
1572
0
    return 0;
1573
282
1574
282
  
if (282
!tab->con[r1].is_row282
)
{0
1575
0
    if (isl_tab_kill_col(tab, tab->con[r1].index) < 0)
1576
0
      return -1;
1577
282
  } else 
if (282
!tab->con[r2].is_row282
)
{0
1578
0
    if (isl_tab_kill_col(tab, tab->con[r2].index) < 0)
1579
0
      return -1;
1580
0
  }
1581
282
1582
282
  
if (282
tab->bmap282
)
{0
1583
0
    tab->bmap = isl_basic_map_add_ineq(tab->bmap, eq);
1584
0
    if (isl_tab_push(tab, isl_tab_undo_bmap_ineq) < 0)
1585
0
      return -1;
1586
0
    isl_seq_neg(eq, eq, 1 + tab->n_var);
1587
0
    tab->bmap = isl_basic_map_add_ineq(tab->bmap, eq);
1588
0
    isl_seq_neg(eq, eq, 1 + tab->n_var);
1589
0
    if (isl_tab_push(tab, isl_tab_undo_bmap_ineq) < 0)
1590
0
      return -1;
1591
0
    
if (0
!tab->bmap0
)
1592
0
      return -1;
1593
0
  }
1594
282
1595
282
  return 0;
1596
282
}
1597
1598
/* Add an inequality to the tableau, resolving violations using
1599
 * restore_lexmin.
1600
 *
1601
 * This function assumes that at least one more row and at least
1602
 * one more element in the constraint array are available in the tableau.
1603
 */
1604
static struct isl_tab *add_lexmin_ineq(struct isl_tab *tab, isl_int *ineq)
1605
10.6k
{
1606
10.6k
  int r;
1607
10.6k
1608
10.6k
  if (!tab)
1609
0
    return NULL;
1610
10.6k
  
if (10.6k
tab->bmap10.6k
)
{0
1611
0
    tab->bmap = isl_basic_map_add_ineq(tab->bmap, ineq);
1612
0
    if (isl_tab_push(tab, isl_tab_undo_bmap_ineq) < 0)
1613
0
      goto error;
1614
0
    
if (0
!tab->bmap0
)
1615
0
      goto error;
1616
0
  }
1617
10.6k
  r = isl_tab_add_row(tab, ineq);
1618
10.6k
  if (r < 0)
1619
0
    goto error;
1620
10.6k
  tab->con[r].is_nonneg = 1;
1621
10.6k
  if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r]) < 0)
1622
0
    goto error;
1623
10.6k
  
if (10.6k
isl_tab_row_is_redundant(tab, tab->con[r].index)10.6k
)
{302
1624
302
    if (isl_tab_mark_redundant(tab, tab->con[r].index) < 0)
1625
0
      goto error;
1626
302
    return tab;
1627
302
  }
1628
10.6k
1629
10.3k
  
if (10.3k
restore_lexmin(tab) < 010.3k
)
1630
0
    goto error;
1631
10.3k
  
if (10.3k
!tab->empty && 10.3k
tab->con[r].is_row9.93k
&&
1632
7.86k
     isl_tab_row_is_redundant(tab, tab->con[r].index))
1633
0
    
if (0
isl_tab_mark_redundant(tab, tab->con[r].index) < 00
)
1634
0
      goto error;
1635
10.3k
  return tab;
1636
0
error:
1637
0
  isl_tab_free(tab);
1638
0
  return NULL;
1639
10.3k
}
1640
1641
/* Check if the coefficients of the parameters are all integral.
1642
 */
1643
static int integer_parameter(struct isl_tab *tab, int row)
1644
10.9k
{
1645
10.9k
  int i;
1646
10.9k
  int col;
1647
10.9k
  unsigned off = 2 + tab->M;
1648
10.9k
1649
41.5k
  for (i = 0; 
i < tab->n_param41.5k
;
++i30.5k
)
{30.9k
1650
30.9k
    /* Eliminated parameter */
1651
30.9k
    if (tab->var[i].is_row)
1652
15.4k
      continue;
1653
15.4k
    col = tab->var[i].index;
1654
15.4k
    if (
!15.4k
isl_int_is_divisible_by15.4k
(tab->mat->row[row][off + col],
1655
15.4k
            tab->mat->row[row][0]))
1656
386
      return 0;
1657
15.4k
  }
1658
15.6k
  
for (i = 0; 10.6k
i < tab->n_div15.6k
;
++i5.03k
)
{5.04k
1659
5.04k
    if (tab->var[tab->n_var - tab->n_div + i].is_row)
1660
1.98k
      continue;
1661
3.06k
    col = tab->var[tab->n_var - tab->n_div + i].index;
1662
3.06k
    if (
!3.06k
isl_int_is_divisible_by3.06k
(tab->mat->row[row][off + col],
1663
3.06k
            tab->mat->row[row][0]))
1664
15
      return 0;
1665
3.06k
  }
1666
10.5k
  return 1;
1667
10.6k
}
1668
1669
/* Check if the coefficients of the non-parameter variables are all integral.
1670
 */
1671
static int integer_variable(struct isl_tab *tab, int row)
1672
802
{
1673
802
  int i;
1674
802
  unsigned off = 2 + tab->M;
1675
802
1676
1.88k
  for (i = tab->n_dead; 
i < tab->n_col1.88k
;
++i1.07k
)
{1.84k
1677
1.84k
    if (tab->col_var[i] >= 0 &&
1678
727
        (tab->col_var[i] < tab->n_param ||
1679
85
         tab->col_var[i] >= tab->n_var - tab->n_div))
1680
699
      continue;
1681
1.14k
    
if (1.14k
!1.14k
isl_int_is_divisible_by1.14k
(tab->mat->row[row][off + i],
1682
1.14k
            tab->mat->row[row][0]))
1683
766
      return 0;
1684
1.14k
  }
1685
36
  return 1;
1686
802
}
1687
1688
/* Check if the constant term is integral.
1689
 */
1690
static int integer_constant(struct isl_tab *tab, int row)
1691
10.9k
{
1692
10.9k
  return isl_int_is_divisible_by(tab->mat->row[row][1],
1693
10.9k
          tab->mat->row[row][0]);
1694
10.9k
}
1695
1696
#define I_CST 1 << 0
1697
#define I_PAR 1 << 1
1698
#define I_VAR 1 << 2
1699
1700
/* Check for next (non-parameter) variable after "var" (first if var == -1)
1701
 * that is non-integer and therefore requires a cut and return
1702
 * the index of the variable.
1703
 * For parametric tableaus, there are three parts in a row,
1704
 * the constant, the coefficients of the parameters and the rest.
1705
 * For each part, we check whether the coefficients in that part
1706
 * are all integral and if so, set the corresponding flag in *f.
1707
 * If the constant and the parameter part are integral, then the
1708
 * current sample value is integral and no cut is required
1709
 * (irrespective of whether the variable part is integral).
1710
 */
1711
static int next_non_integer_var(struct isl_tab *tab, int var, int *f)
1712
4.40k
{
1713
4.40k
  var = var < 0 ? 
tab->n_param4.40k
:
var + 10
;
1714
4.40k
1715
20.4k
  for (; 
var < tab->n_var - tab->n_div20.4k
;
++var16.0k
)
{16.8k
1716
16.8k
    int flags = 0;
1717
16.8k
    int row;
1718
16.8k
    if (!tab->var[var].is_row)
1719
5.83k
      continue;
1720
10.9k
    row = tab->var[var].index;
1721
10.9k
    if (integer_constant(tab, row))
1722
10.3k
      ISL_FL_SET(flags, I_CST);
1723
10.9k
    if (integer_parameter(tab, row))
1724
10.5k
      ISL_FL_SET(flags, I_PAR);
1725
10.9k
    if (
ISL_FL_ISSET10.9k
(flags, I_CST) && 10.9k
ISL_FL_ISSET10.3k
(flags, I_PAR))
1726
10.1k
      continue;
1727
802
    
if (802
integer_variable(tab, row)802
)
1728
36
      ISL_FL_SET(flags, I_VAR);
1729
802
    *f = flags;
1730
802
    return var;
1731
10.9k
  }
1732
3.60k
  return -1;
1733
4.40k
}
1734
1735
/* Check for first (non-parameter) variable that is non-integer and
1736
 * therefore requires a cut and return the corresponding row.
1737
 * For parametric tableaus, there are three parts in a row,
1738
 * the constant, the coefficients of the parameters and the rest.
1739
 * For each part, we check whether the coefficients in that part
1740
 * are all integral and if so, set the corresponding flag in *f.
1741
 * If the constant and the parameter part are integral, then the
1742
 * current sample value is integral and no cut is required
1743
 * (irrespective of whether the variable part is integral).
1744
 */
1745
static int first_non_integer_row(struct isl_tab *tab, int *f)
1746
3.68k
{
1747
3.68k
  int var = next_non_integer_var(tab, -1, f);
1748
3.68k
1749
2.95k
  return var < 0 ? 
-12.95k
:
tab->var[var].index724
;
1750
3.68k
}
1751
1752
/* Add a (non-parametric) cut to cut away the non-integral sample
1753
 * value of the given row.
1754
 *
1755
 * If the row is given by
1756
 *
1757
 *  m r = f + \sum_i a_i y_i
1758
 *
1759
 * then the cut is
1760
 *
1761
 *  c = - {-f/m} + \sum_i {a_i/m} y_i >= 0
1762
 *
1763
 * The big parameter, if any, is ignored, since it is assumed to be big
1764
 * enough to be divisible by any integer.
1765
 * If the tableau is actually a parametric tableau, then this function
1766
 * is only called when all coefficients of the parameters are integral.
1767
 * The cut therefore has zero coefficients for the parameters.
1768
 *
1769
 * The current value is known to be negative, so row_sign, if it
1770
 * exists, is set accordingly.
1771
 *
1772
 * Return the row of the cut or -1.
1773
 */
1774
static int add_cut(struct isl_tab *tab, int row)
1775
401
{
1776
401
  int i;
1777
401
  int r;
1778
401
  isl_int *r_row;
1779
401
  unsigned off = 2 + tab->M;
1780
401
1781
401
  if (isl_tab_extend_cons(tab, 1) < 0)
1782
0
    return -1;
1783
401
  r = isl_tab_allocate_con(tab);
1784
401
  if (r < 0)
1785
0
    return -1;
1786
401
1787
401
  r_row = tab->mat->row[tab->con[r].index];
1788
401
  isl_int_set(r_row[0], tab->mat->row[row][0]);
1789
401
  isl_int_neg(r_row[1], tab->mat->row[row][1]);
1790
401
  isl_int_fdiv_r(r_row[1], r_row[1], tab->mat->row[row][0]);
1791
401
  isl_int_neg(r_row[1], r_row[1]);
1792
401
  if (tab->M)
1793
323
    isl_int_set_si(r_row[2], 0);
1794
4.40k
  for (i = 0; 
i < tab->n_col4.40k
;
++i4.00k
)
1795
4.00k
    isl_int_fdiv_r(r_row[off + i],
1796
401
      tab->mat->row[row][off + i], tab->mat->row[row][0]);
1797
401
1798
401
  tab->con[r].is_nonneg = 1;
1799
401
  if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r]) < 0)
1800
0
    return -1;
1801
401
  
if (401
tab->row_sign401
)
1802
323
    tab->row_sign[tab->con[r].index] = isl_tab_row_neg;
1803
401
1804
401
  return tab->con[r].index;
1805
401
}
1806
1807
0
#define CUT_ALL 1
1808
1.17k
#define CUT_ONE 0
1809
1810
/* Given a non-parametric tableau, add cuts until an integer
1811
 * sample point is obtained or until the tableau is determined
1812
 * to be integer infeasible.
1813
 * As long as there is any non-integer value in the sample point,
1814
 * we add appropriate cuts, if possible, for each of these
1815
 * non-integer values and then resolve the violated
1816
 * cut constraints using restore_lexmin.
1817
 * If one of the corresponding rows is equal to an integral
1818
 * combination of variables/constraints plus a non-integral constant,
1819
 * then there is no way to obtain an integer point and we return
1820
 * a tableau that is marked empty.
1821
 * The parameter cutting_strategy controls the strategy used when adding cuts
1822
 * to remove non-integer points. CUT_ALL adds all possible cuts
1823
 * before continuing the search. CUT_ONE adds only one cut at a time.
1824
 */
1825
static struct isl_tab *cut_to_integer_lexmin(struct isl_tab *tab,
1826
  int cutting_strategy)
1827
1.09k
{
1828
1.09k
  int var;
1829
1.09k
  int row;
1830
1.09k
  int flags;
1831
1.09k
1832
1.09k
  if (!tab)
1833
0
    return NULL;
1834
1.09k
  
if (1.09k
tab->empty1.09k
)
1835
446
    return tab;
1836
1.09k
1837
729
  
while (653
(var = next_non_integer_var(tab, -1, &flags)) != -1729
)
{78
1838
78
    do {
1839
78
      if (
ISL_FL_ISSET78
(flags, I_VAR))
{0
1840
0
        if (isl_tab_mark_empty(tab) < 0)
1841
0
          goto error;
1842
0
        return tab;
1843
0
      }
1844
78
      row = tab->var[var].index;
1845
78
      row = add_cut(tab, row);
1846
78
      if (row < 0)
1847
0
        goto error;
1848
78
      
if (78
cutting_strategy == 78
CUT_ONE78
)
1849
78
        break;
1850
0
    } while ((var = next_non_integer_var(tab, var, &flags)) != -1);
1851
78
    
if (78
restore_lexmin(tab) < 078
)
1852
0
      goto error;
1853
78
    
if (78
tab->empty78
)
1854
2
      break;
1855
78
  }
1856
653
  return tab;
1857
0
error:
1858
0
  isl_tab_free(tab);
1859
0
  return NULL;
1860
653
}
1861
1862
/* Check whether all the currently active samples also satisfy the inequality
1863
 * "ineq" (treated as an equality if eq is set).
1864
 * Remove those samples that do not.
1865
 */
1866
static struct isl_tab *check_samples(struct isl_tab *tab, isl_int *ineq, int eq)
1867
5.87k
{
1868
5.87k
  int i;
1869
5.87k
  isl_int v;
1870
5.87k
1871
5.87k
  if (!tab)
1872
0
    return NULL;
1873
5.87k
1874
5.87k
  
isl_assert5.87k
(tab->mat->ctx, tab->bmap, goto error);5.87k
1875
5.87k
  
isl_assert5.87k
(tab->mat->ctx, tab->samples, goto error);5.87k
1876
5.87k
  
isl_assert5.87k
(tab->mat->ctx, tab->samples->n_col == 1 + tab->n_var, goto error);5.87k
1877
5.87k
1878
5.87k
  
isl_int_init5.87k
(v);5.87k
1879
14.9k
  for (i = tab->n_outside; 
i < tab->n_sample14.9k
;
++i9.05k
)
{9.05k
1880
9.05k
    int sgn;
1881
9.05k
    isl_seq_inner_product(ineq, tab->samples->row[i],
1882
9.05k
          1 + tab->n_var, &v);
1883
9.05k
    sgn = isl_int_sgn(v);
1884
9.05k
    if (
eq ? 9.05k
(sgn == 0)3.60k
:
(sgn >= 0)5.45k
)
1885
6.28k
      continue;
1886
2.76k
    tab = isl_tab_drop_sample(tab, i);
1887
2.76k
    if (!tab)
1888
0
      break;
1889
2.76k
  }
1890
5.87k
  isl_int_clear(v);
1891
5.87k
1892
5.87k
  return tab;
1893
0
error:
1894
0
  isl_tab_free(tab);
1895
0
  return NULL;
1896
5.87k
}
1897
1898
/* Check whether the sample value of the tableau is finite,
1899
 * i.e., either the tableau does not use a big parameter, or
1900
 * all values of the variables are equal to the big parameter plus
1901
 * some constant.  This constant is the actual sample value.
1902
 */
1903
static int sample_is_finite(struct isl_tab *tab)
1904
0
{
1905
0
  int i;
1906
0
1907
0
  if (!tab->M)
1908
0
    return 1;
1909
0
1910
0
  
for (i = 0; 0
i < tab->n_var0
;
++i0
)
{0
1911
0
    int row;
1912
0
    if (!tab->var[i].is_row)
1913
0
      return 0;
1914
0
    row = tab->var[i].index;
1915
0
    if (isl_int_ne(tab->mat->row[row][0], tab->mat->row[row][2]))
1916
0
      return 0;
1917
0
  }
1918
0
  return 1;
1919
0
}
1920
1921
/* Check if the context tableau of sol has any integer points.
1922
 * Leave tab in empty state if no integer point can be found.
1923
 * If an integer point can be found and if moreover it is finite,
1924
 * then it is added to the list of sample values.
1925
 *
1926
 * This function is only called when none of the currently active sample
1927
 * values satisfies the most recently added constraint.
1928
 */
1929
static struct isl_tab *check_integer_feasible(struct isl_tab *tab)
1930
0
{
1931
0
  struct isl_tab_undo *snap;
1932
0
1933
0
  if (!tab)
1934
0
    return NULL;
1935
0
1936
0
  snap = isl_tab_snap(tab);
1937
0
  if (isl_tab_push_basis(tab) < 0)
1938
0
    goto error;
1939
0
1940
0
  
tab = cut_to_integer_lexmin(tab, 0
CUT_ALL0
);
1941
0
  if (!tab)
1942
0
    goto error;
1943
0
1944
0
  
if (0
!tab->empty && 0
sample_is_finite(tab)0
)
{0
1945
0
    struct isl_vec *sample;
1946
0
1947
0
    sample = isl_tab_get_sample_value(tab);
1948
0
1949
0
    if (isl_tab_add_sample(tab, sample) < 0)
1950
0
      goto error;
1951
0
  }
1952
0
1953
0
  
if (0
!tab->empty && 0
isl_tab_rollback(tab, snap) < 00
)
1954
0
    goto error;
1955
0
1956
0
  return tab;
1957
0
error:
1958
0
  isl_tab_free(tab);
1959
0
  return NULL;
1960
0
}
1961
1962
/* Check if any of the currently active sample values satisfies
1963
 * the inequality "ineq" (an equality if eq is set).
1964
 */
1965
static int tab_has_valid_sample(struct isl_tab *tab, isl_int *ineq, int eq)
1966
10.7k
{
1967
10.7k
  int i;
1968
10.7k
  isl_int v;
1969
10.7k
1970
10.7k
  if (!tab)
1971
0
    return -1;
1972
10.7k
1973
10.7k
  
isl_assert10.7k
(tab->mat->ctx, tab->bmap, return -1);10.7k
1974
10.7k
  
isl_assert10.7k
(tab->mat->ctx, tab->samples, return -1);10.7k
1975
10.7k
  
isl_assert10.7k
(tab->mat->ctx, tab->samples->n_col == 1 + tab->n_var, return -1);10.7k
1976
10.7k
1977
10.7k
  
isl_int_init10.7k
(v);10.7k
1978
17.8k
  for (i = tab->n_outside; 
i < tab->n_sample17.8k
;
++i7.12k
)
{10.7k
1979
10.7k
    int sgn;
1980
10.7k
    isl_seq_inner_product(ineq, tab->samples->row[i],
1981
10.7k
          1 + tab->n_var, &v);
1982
10.7k
    sgn = isl_int_sgn(v);
1983
10.7k
    if (
eq ? 10.7k
(sgn == 0)3.54k
:
(sgn >= 0)7.16k
)
1984
3.58k
      break;
1985
10.7k
  }
1986
10.7k
  isl_int_clear(v);
1987
10.7k
1988
10.7k
  return i < tab->n_sample;
1989
10.7k
}
1990
1991
/* Insert a div specified by "div" to the tableau "tab" at position "pos" and
1992
 * return isl_bool_true if the div is obviously non-negative.
1993
 */
1994
static isl_bool context_tab_insert_div(struct isl_tab *tab, int pos,
1995
  __isl_keep isl_vec *div,
1996
  isl_stat (*add_ineq)(void *user, isl_int *), void *user)
1997
386
{
1998
386
  int i;
1999
386
  int r;
2000
386
  struct isl_mat *samples;
2001
386
  int nonneg;
2002
386
2003
386
  r = isl_tab_insert_div(tab, pos, div, add_ineq, user);
2004
386
  if (r < 0)
2005
0
    return isl_bool_error;
2006
386
  nonneg = tab->var[r].is_nonneg;
2007
386
  tab->var[r].frozen = 1;
2008
386
2009
386
  samples = isl_mat_extend(tab->samples,
2010
386
      tab->n_sample, 1 + tab->n_var);
2011
386
  tab->samples = samples;
2012
386
  if (!samples)
2013
0
    return isl_bool_error;
2014
1.57k
  
for (i = tab->n_outside; 386
i < samples->n_row1.57k
;
++i1.18k
)
{1.18k
2015
1.18k
    isl_seq_inner_product(div->el + 1, samples->row[i],
2016
1.18k
      div->size - 1, &samples->row[i][samples->n_col - 1]);
2017
1.18k
    isl_int_fdiv_q(samples->row[i][samples->n_col - 1],
2018
1.18k
             samples->row[i][samples->n_col - 1], div->el[0]);
2019
1.18k
  }
2020
386
  tab->samples = isl_mat_move_cols(tab->samples, 1 + pos,
2021
386
          1 + tab->n_var - 1, 1);
2022
386
  if (!tab->samples)
2023
0
    return isl_bool_error;
2024
386
2025
386
  return nonneg;
2026
386
}
2027
2028
/* Add a div specified by "div" to both the main tableau and
2029
 * the context tableau.  In case of the main tableau, we only
2030
 * need to add an extra div.  In the context tableau, we also
2031
 * need to express the meaning of the div.
2032
 * Return the index of the div or -1 if anything went wrong.
2033
 *
2034
 * The new integer division is added before any unknown integer
2035
 * divisions in the context to ensure that it does not get
2036
 * equated to some linear combination involving unknown integer
2037
 * divisions.
2038
 */
2039
static int add_div(struct isl_tab *tab, struct isl_context *context,
2040
  __isl_keep isl_vec *div)
2041
386
{
2042
386
  int r;
2043
386
  int pos;
2044
386
  isl_bool nonneg;
2045
386
  struct isl_tab *context_tab = context->op->peek_tab(context);
2046
386
2047
386
  if (
!tab || 386
!context_tab386
)
2048
0
    goto error;
2049
386
2050
386
  pos = context_tab->n_var - context->n_unknown;
2051
386
  if ((nonneg = context->op->insert_div(context, pos, div)) < 0)
2052
0
    goto error;
2053
386
2054
386
  
if (386
!context->op->is_ok(context)386
)
2055
0
    goto error;
2056
386
2057
386
  pos = tab->n_var - context->n_unknown;
2058
386
  if (isl_tab_extend_vars(tab, 1) < 0)
2059
0
    goto error;
2060
386
  r = isl_tab_insert_var(tab, pos);
2061
386
  if (r < 0)
2062
0
    goto error;
2063
386
  
if (386
nonneg386
)
2064
47
    tab->var[r].is_nonneg = 1;
2065
386
  tab->var[r].frozen = 1;
2066
386
  tab->n_div++;
2067
386
2068
386
  return tab->n_div - 1 - context->n_unknown;
2069
0
error:
2070
0
  context->op->invalidate(context);
2071
0
  return -1;
2072
386
}
2073
2074
static int find_div(struct isl_tab *tab, isl_int *div, isl_int denom)
2075
401
{
2076
401
  int i;
2077
401
  unsigned total = isl_basic_map_total_dim(tab->bmap);
2078
401
2079
1.56k
  for (i = 0; 
i < tab->bmap->n_div1.56k
;
++i1.16k
)
{1.17k
2080
1.17k
    if (isl_int_ne(tab->bmap->div[i][0], denom))
2081
847
      continue;
2082
329
    
if (329
!isl_seq_eq(tab->bmap->div[i] + 1, div, 1 + total)329
)
2083
314
      continue;
2084
15
    return i;
2085
329
  }
2086
386
  return -1;
2087
401
}
2088
2089
/* Return the index of a div that corresponds to "div".
2090
 * We first check if we already have such a div and if not, we create one.
2091
 */
2092
static int get_div(struct isl_tab *tab, struct isl_context *context,
2093
  struct isl_vec *div)
2094
401
{
2095
401
  int d;
2096
401
  struct isl_tab *context_tab = context->op->peek_tab(context);
2097
401
2098
401
  if (!context_tab)
2099
0
    return -1;
2100
401
2101
401
  d = find_div(context_tab, div->el + 1, div->el[0]);
2102
401
  if (d != -1)
2103
15
    return d;
2104
401
2105
386
  return add_div(tab, context, div);
2106
401
}
2107
2108
/* Add a parametric cut to cut away the non-integral sample value
2109
 * of the give row.
2110
 * Let a_i be the coefficients of the constant term and the parameters
2111
 * and let b_i be the coefficients of the variables or constraints
2112
 * in basis of the tableau.
2113
 * Let q be the div q = floor(\sum_i {-a_i} y_i).
2114
 *
2115
 * The cut is expressed as
2116
 *
2117
 *  c = \sum_i -{-a_i} y_i + \sum_i {b_i} x_i + q >= 0
2118
 *
2119
 * If q did not already exist in the context tableau, then it is added first.
2120
 * If q is in a column of the main tableau then the "+ q" can be accomplished
2121
 * by setting the corresponding entry to the denominator of the constraint.
2122
 * If q happens to be in a row of the main tableau, then the corresponding
2123
 * row needs to be added instead (taking care of the denominators).
2124
 * Note that this is very unlikely, but perhaps not entirely impossible.
2125
 *
2126
 * The current value of the cut is known to be negative (or at least
2127
 * non-positive), so row_sign is set accordingly.
2128
 *
2129
 * Return the row of the cut or -1.
2130
 */
2131
static int add_parametric_cut(struct isl_tab *tab, int row,
2132
  struct isl_context *context)
2133
365
{
2134
365
  struct isl_vec *div;
2135
365
  int d;
2136
365
  int i;
2137
365
  int r;
2138
365
  isl_int *r_row;
2139
365
  int col;
2140
365
  int n;
2141
365
  unsigned off = 2 + tab->M;
2142
365
2143
365
  if (!context)
2144
0
    return -1;
2145
365
2146
365
  div = get_row_parameter_div(tab, row);
2147
365
  if (!div)
2148
0
    return -1;
2149
365
2150
365
  n = tab->n_div - context->n_unknown;
2151
365
  d = context->op->get_div(context, tab, div);
2152
365
  isl_vec_free(div);
2153
365
  if (d < 0)
2154
0
    return -1;
2155
365
2156
365
  
if (365
isl_tab_extend_cons(tab, 1) < 0365
)
2157
0
    return -1;
2158
365
  r = isl_tab_allocate_con(tab);
2159
365
  if (r < 0)
2160
0
    return -1;
2161
365
2162
365
  r_row = tab->mat->row[tab->con[r].index];
2163
365
  isl_int_set(r_row[0], tab->mat->row[row][0]);
2164
365
  isl_int_neg(r_row[1], tab->mat->row[row][1]);
2165
365
  isl_int_fdiv_r(r_row[1], r_row[1], tab->mat->row[row][0]);
2166
365
  isl_int_neg(r_row[1], r_row[1]);
2167
365
  if (tab->M)
2168
365
    isl_int_set_si(r_row[2], 0);
2169
1.13k
  for (i = 0; 
i < tab->n_param1.13k
;
++i773
)
{773
2170
773
    if (tab->var[i].is_row)
2171
14
      continue;
2172
759
    col = tab->var[i].index;
2173
759
    isl_int_neg(r_row[off + col], tab->mat->row[row][off + col]);
2174
759
    isl_int_fdiv_r(r_row[off + col], r_row[off + col],
2175
759
        tab->mat->row[row][0]);
2176
759
    isl_int_neg(r_row[off + col], r_row[off + col]);
2177
759
  }
2178
1.86k
  for (i = 0; 
i < tab->n_div1.86k
;
++i1.49k
)
{1.49k
2179
1.49k
    if (tab->var[tab->n_var - tab->n_div + i].is_row)
2180
271
      continue;
2181
1.22k
    col = tab->var[tab->n_var - tab->n_div + i].index;
2182
1.22k
    isl_int_neg(r_row[off + col], tab->mat->row[row][off + col]);
2183
1.22k
    isl_int_fdiv_r(r_row[off + col], r_row[off + col],
2184
1.22k
        tab->mat->row[row][0]);
2185
1.22k
    isl_int_neg(r_row[off + col], r_row[off + col]);
2186
1.22k
  }
2187
3.76k
  for (i = 0; 
i < tab->n_col3.76k
;
++i3.39k
)
{3.39k
2188
3.39k
    if (tab->col_var[i] >= 0 &&
2189
1.98k
        (tab->col_var[i] < tab->n_param ||
2190
1.22k
         tab->col_var[i] >= tab->n_var - tab->n_div))
2191
1.98k
      continue;
2192
1.41k
    
isl_int_fdiv_r1.41k
(r_row[off + i],1.41k
2193
1.41k
      tab->mat->row[row][off + i], tab->mat->row[row][0]);
2194
1.41k
  }
2195
365
  if (
tab->var[tab->n_var - tab->n_div + d].is_row365
)
{1
2196
1
    isl_int gcd;
2197
1
    int d_row = tab->var[tab->n_var - tab->n_div + d].index;
2198
1
    isl_int_init(gcd);
2199
1
    isl_int_gcd(gcd, tab->mat->row[d_row][0], r_row[0]);
2200
1
    isl_int_divexact(r_row[0], r_row[0], gcd);
2201
1
    isl_int_divexact(gcd, tab->mat->row[d_row][0], gcd);
2202
1
    isl_seq_combine(r_row + 1, gcd, r_row + 1,
2203
1
        r_row[0], tab->mat->row[d_row] + 1,
2204
1
        off - 1 + tab->n_col);
2205
1
    isl_int_mul(r_row[0], r_row[0], tab->mat->row[d_row][0]);
2206
1
    isl_int_clear(gcd);
2207
364
  } else {
2208
364
    col = tab->var[tab->n_var - tab->n_div + d].index;
2209
364
    isl_int_set(r_row[off + col], tab->mat->row[row][0]);
2210
364
  }
2211
365
2212
365
  tab->con[r].is_nonneg = 1;
2213
365
  if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r]) < 0)
2214
0
    return -1;
2215
365
  
if (365
tab->row_sign365
)
2216
365
    tab->row_sign[tab->con[r].index] = isl_tab_row_neg;
2217
365
2218
365
  row = tab->con[r].index;
2219
365
2220
365
  if (
d >= n && 365
context->op->detect_equalities(context, tab) < 0363
)
2221
0
    return -1;
2222
365
2223
365
  return row;
2224
365
}
2225
2226
/* Construct a tableau for bmap that can be used for computing
2227
 * the lexicographic minimum (or maximum) of bmap.
2228
 * If not NULL, then dom is the domain where the minimum
2229
 * should be computed.  In this case, we set up a parametric
2230
 * tableau with row signs (initialized to "unknown").
2231
 * If M is set, then the tableau will use a big parameter.
2232
 * If max is set, then a maximum should be computed instead of a minimum.
2233
 * This means that for each variable x, the tableau will contain the variable
2234
 * x' = M - x, rather than x' = M + x.  This in turn means that the coefficient
2235
 * of the variables in all constraints are negated prior to adding them
2236
 * to the tableau.
2237
 */
2238
static struct isl_tab *tab_for_lexmin(struct isl_basic_map *bmap,
2239
  struct isl_basic_set *dom, unsigned M, int max)
2240
3.20k
{
2241
3.20k
  int i;
2242
3.20k
  struct isl_tab *tab;
2243
3.20k
  unsigned n_var;
2244
3.20k
  unsigned o_var;
2245
3.20k
2246
3.20k
  tab = isl_tab_alloc(bmap->ctx, 2 * bmap->n_eq + bmap->n_ineq + 1,
2247
3.20k
          isl_basic_map_total_dim(bmap), M);
2248
3.20k
  if (!tab)
2249
0
    return NULL;
2250
3.20k
2251
3.20k
  
tab->rational = 3.20k
ISL_F_ISSET3.20k
(bmap, ISL_BASIC_MAP_RATIONAL);
2252
3.20k
  if (
dom3.20k
)
{2.80k
2253
2.80k
    tab->n_param = isl_basic_set_total_dim(dom) - dom->n_div;
2254
2.80k
    tab->n_div = dom->n_div;
2255
2.80k
    tab->row_sign = isl_calloc_array(bmap->ctx,
2256
2.80k
          enum isl_tab_row_sign, tab->mat->n_row);
2257
2.80k
    if (
tab->mat->n_row && 2.80k
!tab->row_sign2.80k
)
2258
0
      goto error;
2259
2.80k
  }
2260
3.20k
  
if (3.20k
ISL_F_ISSET3.20k
(bmap, ISL_BASIC_MAP_EMPTY))
{0
2261
0
    if (isl_tab_mark_empty(tab) < 0)
2262
0
      goto error;
2263
0
    return tab;
2264
0
  }
2265
3.20k
2266
14.4k
  
for (i = tab->n_param; 3.20k
i < tab->n_var - tab->n_div14.4k
;
++i11.2k
)
{11.2k
2267
11.2k
    tab->var[i].is_nonneg = 1;
2268
11.2k
    tab->var[i].frozen = 1;
2269
11.2k
  }
2270
3.20k
  o_var = 1 + tab->n_param;
2271
3.20k
  n_var = tab->n_var - tab->n_param - tab->n_div;
2272
12.6k
  for (i = 0; 
i < bmap->n_eq12.6k
;
++i9.47k
)
{9.47k
2273
9.47k
    if (max)
2274
6.76k
      isl_seq_neg(bmap->eq[i] + o_var,
2275
6.76k
            bmap->eq[i] + o_var, n_var);
2276
9.47k
    tab = add_lexmin_valid_eq(tab, bmap->eq[i]);
2277
9.47k
    if (max)
2278
6.76k
      isl_seq_neg(bmap->eq[i] + o_var,
2279
6.76k
            bmap->eq[i] + o_var, n_var);
2280
9.47k
    if (
!tab || 9.47k
tab->empty9.47k
)
2281
0
      return tab;
2282
9.47k
  }
2283
3.20k
  
if (3.20k
bmap->n_eq && 3.20k
restore_lexmin(tab) < 02.53k
)
2284
0
    goto error;
2285
13.1k
  
for (i = 0; 3.20k
i < bmap->n_ineq13.1k
;
++i9.94k
)
{9.94k
2286
9.94k
    if (max)
2287
4.65k
      isl_seq_neg(bmap->ineq[i] + o_var,
2288
4.65k
            bmap->ineq[i] + o_var, n_var);
2289
9.94k
    tab = add_lexmin_ineq(tab, bmap->ineq[i]);
2290
9.94k
    if (max)
2291
4.65k
      isl_seq_neg(bmap->ineq[i] + o_var,
2292
4.65k
            bmap->ineq[i] + o_var, n_var);
2293
9.94k
    if (
!tab || 9.94k
tab->empty9.94k
)
2294
0
      return tab;
2295
9.94k
  }
2296
3.20k
  return tab;
2297
0
error:
2298
0
  isl_tab_free(tab);
2299
0
  return NULL;
2300
3.20k
}
2301
2302
/* Given a main tableau where more than one row requires a split,
2303
 * determine and return the "best" row to split on.
2304
 *
2305
 * Given two rows in the main tableau, if the inequality corresponding
2306
 * to the first row is redundant with respect to that of the second row
2307
 * in the current tableau, then it is better to split on the second row,
2308
 * since in the positive part, both rows will be positive.
2309
 * (In the negative part a pivot will have to be performed and just about
2310
 * anything can happen to the sign of the other row.)
2311
 *
2312
 * As a simple heuristic, we therefore select the row that makes the most
2313
 * of the other rows redundant.
2314
 *
2315
 * Perhaps it would also be useful to look at the number of constraints
2316
 * that conflict with any given constraint.
2317
 *
2318
 * best is the best row so far (-1 when we have not found any row yet).
2319
 * best_r is the number of other rows made redundant by row best.
2320
 * When best is still -1, bset_r is meaningless, but it is initialized
2321
 * to some arbitrary value (0) anyway.  Without this redundant initialization
2322
 * valgrind may warn about uninitialized memory accesses when isl
2323
 * is compiled with some versions of gcc.
2324
 */
2325
static int best_split(struct isl_tab *tab, struct isl_tab *context_tab)
2326
235
{
2327
235
  struct isl_tab_undo *snap;
2328
235
  int split;
2329
235
  int row;
2330
235
  int best = -1;
2331
235
  int best_r = 0;
2332
235
2333
235
  if (isl_tab_extend_cons(context_tab, 2) < 0)
2334
0
    return -1;
2335
235
2336
235
  snap = isl_tab_snap(context_tab);
2337
235
2338
2.06k
  for (split = tab->n_redundant; 
split < tab->n_row2.06k
;
++split1.82k
)
{1.82k
2339
1.82k
    struct isl_tab_undo *snap2;
2340
1.82k
    struct isl_vec *ineq = NULL;
2341
1.82k
    int r = 0;
2342
1.82k
    int ok;
2343
1.82k
2344
1.82k
    if (!isl_tab_var_from_row(tab, split)->is_nonneg)
2345
34
      continue;
2346
1.79k
    
if (1.79k
tab->row_sign[split] != isl_tab_row_any1.79k
)
2347
1.22k
      continue;
2348
1.79k
2349
571
    ineq = get_row_parameter_ineq(tab, split);
2350
571
    if (!ineq)
2351
0
      return -1;
2352
571
    ok = isl_tab_add_ineq(context_tab, ineq->el) >= 0;
2353
571
    isl_vec_free(ineq);
2354
571
    if (!ok)
2355
0
      return -1;
2356
571
2357
571
    snap2 = isl_tab_snap(context_tab);
2358
571
2359
5.23k
    for (row = tab->n_redundant; 
row < tab->n_row5.23k
;
++row4.66k
)
{4.66k
2360
4.66k
      struct isl_tab_var *var;
2361
4.66k
2362
4.66k
      if (row == split)
2363
571
        continue;
2364
4.09k
      
if (4.09k
!isl_tab_var_from_row(tab, row)->is_nonneg4.09k
)
2365
74
        continue;
2366
4.02k
      
if (4.02k
tab->row_sign[row] != isl_tab_row_any4.02k
)
2367
2.99k
        continue;
2368
4.02k
2369
1.02k
      ineq = get_row_parameter_ineq(tab, row);
2370
1.02k
      if (!ineq)
2371
0
        return -1;
2372
1.02k
      ok = isl_tab_add_ineq(context_tab, ineq->el) >= 0;
2373
1.02k
      isl_vec_free(ineq);
2374
1.02k
      if (!ok)
2375
0
        return -1;
2376
1.02k
      var = &context_tab->con[context_tab->n_con - 1];
2377
1.02k
      if (!context_tab->empty &&
2378
1.02k
          !isl_tab_min_at_most_neg_one(context_tab, var))
2379
47
        r++;
2380
1.02k
      if (isl_tab_rollback(context_tab, snap2) < 0)
2381
0
        return -1;
2382
1.02k
    }
2383
571
    
if (571
best == -1 || 571
r > best_r336
)
{253
2384
253
      best = split;
2385
253
      best_r = r;
2386
253
    }
2387
571
    if (isl_tab_rollback(context_tab, snap) < 0)
2388
0
      return -1;
2389
571
  }
2390
235
2391
235
  return best;
2392
235
}
2393
2394
static struct isl_basic_set *context_lex_peek_basic_set(
2395
  struct isl_context *context)
2396
0
{
2397
0
  struct isl_context_lex *clex = (struct isl_context_lex *)context;
2398
0
  if (!clex->tab)
2399
0
    return NULL;
2400
0
  return isl_tab_peek_bset(clex->tab);
2401
0
}
2402
2403
static struct isl_tab *context_lex_peek_tab(struct isl_context *context)
2404
0
{
2405
0
  struct isl_context_lex *clex = (struct isl_context_lex *)context;
2406
0
  return clex->tab;
2407
0
}
2408
2409
static void context_lex_add_eq(struct isl_context *context, isl_int *eq,
2410
    int check, int update)
2411
0
{
2412
0
  struct isl_context_lex *clex = (struct isl_context_lex *)context;
2413
0
  if (isl_tab_extend_cons(clex->tab, 2) < 0)
2414
0
    goto error;
2415
0
  
if (0
add_lexmin_eq(clex->tab, eq) < 00
)
2416
0
    goto error;
2417
0
  
if (0
check0
)
{0
2418
0
    int v = tab_has_valid_sample(clex->tab, eq, 1);
2419
0
    if (v < 0)
2420
0
      goto error;
2421
0
    
if (0
!v0
)
2422
0
      clex->tab = check_integer_feasible(clex->tab);
2423
0
  }
2424
0
  
if (0
update0
)
2425
0
    clex->tab = check_samples(clex->tab, eq, 1);
2426
0
  return;
2427
0
error:
2428
0
  isl_tab_free(clex->tab);
2429
0
  clex->tab = NULL;
2430
0
}
2431
2432
static void context_lex_add_ineq(struct isl_context *context, isl_int *ineq,
2433
    int check, int update)
2434
0
{
2435
0
  struct isl_context_lex *clex = (struct isl_context_lex *)context;
2436
0
  if (isl_tab_extend_cons(clex->tab, 1) < 0)
2437
0
    goto error;
2438
0
  clex->tab = add_lexmin_ineq(clex->tab, ineq);
2439
0
  if (
check0
)
{0
2440
0
    int v = tab_has_valid_sample(clex->tab, ineq, 0);
2441
0
    if (v < 0)
2442
0
      goto error;
2443
0
    
if (0
!v0
)
2444
0
      clex->tab = check_integer_feasible(clex->tab);
2445
0
  }
2446
0
  
if (0
update0
)
2447
0
    clex->tab = check_samples(clex->tab, ineq, 0);
2448
0
  return;
2449
0
error:
2450
0
  isl_tab_free(clex->tab);
2451
0
  clex->tab = NULL;
2452
0
}
2453
2454
static isl_stat context_lex_add_ineq_wrap(void *user, isl_int *ineq)
2455
0
{
2456
0
  struct isl_context *context = (struct isl_context *)user;
2457
0
  context_lex_add_ineq(context, ineq, 0, 0);
2458
0
  return context->op->is_ok(context) ? 
isl_stat_ok0
:
isl_stat_error0
;
2459
0
}
2460
2461
/* Check which signs can be obtained by "ineq" on all the currently
2462
 * active sample values.  See row_sign for more information.
2463
 */
2464
static enum isl_tab_row_sign tab_ineq_sign(struct isl_tab *tab, isl_int *ineq,
2465
  int strict)
2466
5.43k
{
2467
5.43k
  int i;
2468
5.43k
  int sgn;
2469
5.43k
  isl_int tmp;
2470
5.43k
  enum isl_tab_row_sign res = isl_tab_row_unknown;
2471
5.43k
2472
5.43k
  isl_assert(tab->mat->ctx, tab->samples, return isl_tab_row_unknown);
2473
5.43k
  
isl_assert5.43k
(tab->mat->ctx, tab->samples->n_col == 1 + tab->n_var,5.43k
2474
5.43k
      return isl_tab_row_unknown);
2475
5.43k
2476
5.43k
  
isl_int_init5.43k
(tmp);5.43k
2477
14.3k
  for (i = tab->n_outside; 
i < tab->n_sample14.3k
;
++i8.92k
)
{9.31k
2478
9.31k
    isl_seq_inner_product(tab->samples->row[i], ineq,
2479
9.31k
          1 + tab->n_var, &tmp);
2480
9.31k
    sgn = isl_int_sgn(tmp);
2481
9.31k
    if (
sgn > 0 || 9.31k
(sgn == 0 && 3.87k
strict2.23k
))
{7.46k
2482
7.46k
      if (res == isl_tab_row_unknown)
2483
4.27k
        res = isl_tab_row_pos;
2484
7.46k
      if (res == isl_tab_row_neg)
2485
194
        res = isl_tab_row_any;
2486
7.46k
    }
2487
9.31k
    if (
sgn < 09.31k
)
{1.64k
2488
1.64k
      if (res == isl_tab_row_unknown)
2489
1.12k
        res = isl_tab_row_neg;
2490
1.64k
      if (res == isl_tab_row_pos)
2491
195
        res = isl_tab_row_any;
2492
1.64k
    }
2493
9.31k
    if (res == isl_tab_row_any)
2494
389
      break;
2495
9.31k
  }
2496
5.43k
  isl_int_clear(tmp);
2497
5.43k
2498
5.43k
  return res;
2499
5.43k
}
2500
2501
static enum isl_tab_row_sign context_lex_ineq_sign(struct isl_context *context,
2502
      isl_int *ineq, int strict)
2503
0
{
2504
0
  struct isl_context_lex *clex = (struct isl_context_lex *)context;
2505
0
  return tab_ineq_sign(clex->tab, ineq, strict);
2506
0
}
2507
2508
/* Check whether "ineq" can be added to the tableau without rendering
2509
 * it infeasible.
2510
 */
2511
static int context_lex_test_ineq(struct isl_context *context, isl_int *ineq)
2512
0
{
2513
0
  struct isl_context_lex *clex = (struct isl_context_lex *)context;
2514
0
  struct isl_tab_undo *snap;
2515
0
  int feasible;
2516
0
2517
0
  if (!clex->tab)
2518
0
    return -1;
2519
0
2520
0
  
if (0
isl_tab_extend_cons(clex->tab, 1) < 00
)
2521
0
    return -1;
2522
0
2523
0
  snap = isl_tab_snap(clex->tab);
2524
0
  if (isl_tab_push_basis(clex->tab) < 0)
2525
0
    return -1;
2526
0
  clex->tab = add_lexmin_ineq(clex->tab, ineq);
2527
0
  clex->tab = check_integer_feasible(clex->tab);
2528
0
  if (!clex->tab)
2529
0
    return -1;
2530
0
  feasible = !clex->tab->empty;
2531
0
  if (isl_tab_rollback(clex->tab, snap) < 0)
2532
0
    return -1;
2533
0
2534
0
  return feasible;
2535
0
}
2536
2537
static int context_lex_get_div(struct isl_context *context, struct isl_tab *tab,
2538
    struct isl_vec *div)
2539
0
{
2540
0
  return get_div(tab, context, div);
2541
0
}
2542
2543
/* Insert a div specified by "div" to the context tableau at position "pos" and
2544
 * return isl_bool_true if the div is obviously non-negative.
2545
 * context_tab_add_div will always return isl_bool_true, because all variables
2546
 * in a isl_context_lex tableau are non-negative.
2547
 * However, if we are using a big parameter in the context, then this only
2548
 * reflects the non-negativity of the variable used to _encode_ the
2549
 * div, i.e., div' = M + div, so we can't draw any conclusions.
2550
 */
2551
static isl_bool context_lex_insert_div(struct isl_context *context, int pos,
2552
  __isl_keep isl_vec *div)
2553
0
{
2554
0
  struct isl_context_lex *clex = (struct isl_context_lex *)context;
2555
0
  isl_bool nonneg;
2556
0
  nonneg = context_tab_insert_div(clex->tab, pos, div,
2557
0
          context_lex_add_ineq_wrap, context);
2558
0
  if (nonneg < 0)
2559
0
    return isl_bool_error;
2560
0
  
if (0
clex->tab->M0
)
2561
0
    return isl_bool_false;
2562
0
  return nonneg;
2563
0
}
2564
2565
static int context_lex_detect_equalities(struct isl_context *context,
2566
    struct isl_tab *tab)
2567
0
{
2568
0
  return 0;
2569
0
}
2570
2571
static int context_lex_best_split(struct isl_context *context,
2572
    struct isl_tab *tab)
2573
0
{
2574
0
  struct isl_context_lex *clex = (struct isl_context_lex *)context;
2575
0
  struct isl_tab_undo *snap;
2576
0
  int r;
2577
0
2578
0
  snap = isl_tab_snap(clex->tab);
2579
0
  if (isl_tab_push_basis(clex->tab) < 0)
2580
0
    return -1;
2581
0
  r = best_split(tab, clex->tab);
2582
0
2583
0
  if (
r >= 0 && 0
isl_tab_rollback(clex->tab, snap) < 00
)
2584
0
    return -1;
2585
0
2586
0
  return r;
2587
0
}
2588
2589
static int context_lex_is_empty(struct isl_context *context)
2590
0
{
2591
0
  struct isl_context_lex *clex = (struct isl_context_lex *)context;
2592
0
  if (!clex->tab)
2593
0
    return -1;
2594
0
  return clex->tab->empty;
2595
0
}
2596
2597
static void *context_lex_save(struct isl_context *context)
2598
0
{
2599
0
  struct isl_context_lex *clex = (struct isl_context_lex *)context;
2600
0
  struct isl_tab_undo *snap;
2601
0
2602
0
  snap = isl_tab_snap(clex->tab);
2603
0
  if (isl_tab_push_basis(clex->tab) < 0)
2604
0
    return NULL;
2605
0
  
if (0
isl_tab_save_samples(clex->tab) < 00
)
2606
0
    return NULL;
2607
0
2608
0
  return snap;
2609
0
}
2610
2611
static void context_lex_restore(struct isl_context *context, void *save)
2612
0
{
2613
0
  struct isl_context_lex *clex = (struct isl_context_lex *)context;
2614
0
  if (
isl_tab_rollback(clex->tab, (struct isl_tab_undo *)save) < 00
)
{0
2615
0
    isl_tab_free(clex->tab);
2616
0
    clex->tab = NULL;
2617
0
  }
2618
0
}
2619
2620
static void context_lex_discard(void *save)
2621
0
{
2622
0
}
2623
2624
static int context_lex_is_ok(struct isl_context *context)
2625
0
{
2626
0
  struct isl_context_lex *clex = (struct isl_context_lex *)context;
2627
0
  return !!clex->tab;
2628
0
}
2629
2630
/* For each variable in the context tableau, check if the variable can
2631
 * only attain non-negative values.  If so, mark the parameter as non-negative
2632
 * in the main tableau.  This allows for a more direct identification of some
2633
 * cases of violated constraints.
2634
 */
2635
static struct isl_tab *tab_detect_nonnegative_parameters(struct isl_tab *tab,
2636
  struct isl_tab *context_tab)
2637
2.80k
{
2638
2.80k
  int i;
2639
2.80k
  struct isl_tab_undo *snap;
2640
2.80k
  struct isl_vec *ineq = NULL;
2641
2.80k
  struct isl_tab_var *var;
2642
2.80k
  int n;
2643
2.80k
2644
2.80k
  if (context_tab->n_var == 0)
2645
585
    return tab;
2646
2.80k
2647
2.22k
  ineq = isl_vec_alloc(tab->mat->ctx, 1 + context_tab->n_var);
2648
2.22k
  if (!ineq)
2649
0
    goto error;
2650
2.22k
2651
2.22k
  
if (2.22k
isl_tab_extend_cons(context_tab, 1) < 02.22k
)
2652
0
    goto error;
2653
2.22k
2654
2.22k
  snap = isl_tab_snap(context_tab);
2655
2.22k
2656
2.22k
  n = 0;
2657
2.22k
  isl_seq_clr(ineq->el, ineq->size);
2658
10.5k
  for (i = 0; 
i < context_tab->n_var10.5k
;
++i8.28k
)
{8.28k
2659
8.28k
    isl_int_set_si(ineq->el[1 + i], 1);
2660
8.28k
    if (isl_tab_add_ineq(context_tab, ineq->el) < 0)
2661
0
      goto error;
2662
8.28k
    var = &context_tab->con[context_tab->n_con - 1];
2663
8.28k
    if (!context_tab->empty &&
2664
8.16k
        
!isl_tab_min_at_most_neg_one(context_tab, var)8.16k
)
{5.58k
2665
5.58k
      int j = i;
2666
5.58k
      if (i >= tab->n_param)
2667
21
        j = i - tab->n_param + tab->n_var - tab->n_div;
2668
5.58k
      tab->var[j].is_nonneg = 1;
2669
5.58k
      n++;
2670
5.58k
    }
2671
8.28k
    isl_int_set_si(ineq->el[1 + i], 0);
2672
8.28k
    if (isl_tab_rollback(context_tab, snap) < 0)
2673
0
      goto error;
2674
8.28k
  }
2675
2.22k
2676
2.22k
  
if (2.22k
context_tab->M && 2.22k
n == context_tab->n_var0
)
{0
2677
0
    context_tab->mat = isl_mat_drop_cols(context_tab->mat, 2, 1);
2678
0
    context_tab->M = 0;
2679
0
  }
2680
2.22k
2681
2.22k
  isl_vec_free(ineq);
2682
2.22k
  return tab;
2683
0
error:
2684
0
  isl_vec_free(ineq);
2685
0
  isl_tab_free(tab);
2686
0
  return NULL;
2687
2.22k
}
2688
2689
static struct isl_tab *context_lex_detect_nonnegative_parameters(
2690
  struct isl_context *context, struct isl_tab *tab)
2691
0
{
2692
0
  struct isl_context_lex *clex = (struct isl_context_lex *)context;
2693
0
  struct isl_tab_undo *snap;
2694
0
2695
0
  if (!tab)
2696
0
    return NULL;
2697
0
2698
0
  snap = isl_tab_snap(clex->tab);
2699
0
  if (isl_tab_push_basis(clex->tab) < 0)
2700
0
    goto error;
2701
0
2702
0
  tab = tab_detect_nonnegative_parameters(tab, clex->tab);
2703
0
2704
0
  if (isl_tab_rollback(clex->tab, snap) < 0)
2705
0
    goto error;
2706
0
2707
0
  return tab;
2708
0
error:
2709
0
  isl_tab_free(tab);
2710
0
  return NULL;
2711
0
}
2712
2713
static void context_lex_invalidate(struct isl_context *context)
2714
0
{
2715
0
  struct isl_context_lex *clex = (struct isl_context_lex *)context;
2716
0
  isl_tab_free(clex->tab);
2717
0
  clex->tab = NULL;
2718
0
}
2719
2720
static __isl_null struct isl_context *context_lex_free(
2721
  struct isl_context *context)
2722
0
{
2723
0
  struct isl_context_lex *clex = (struct isl_context_lex *)context;
2724
0
  isl_tab_free(clex->tab);
2725
0
  free(clex);
2726
0
2727
0
  return NULL;
2728
0
}
2729
2730
struct isl_context_op isl_context_lex_op = {
2731
  context_lex_detect_nonnegative_parameters,
2732
  context_lex_peek_basic_set,
2733
  context_lex_peek_tab,
2734
  context_lex_add_eq,
2735
  context_lex_add_ineq,
2736
  context_lex_ineq_sign,
2737
  context_lex_test_ineq,
2738
  context_lex_get_div,
2739
  context_lex_insert_div,
2740
  context_lex_detect_equalities,
2741
  context_lex_best_split,
2742
  context_lex_is_empty,
2743
  context_lex_is_ok,
2744
  context_lex_save,
2745
  context_lex_restore,
2746
  context_lex_discard,
2747
  context_lex_invalidate,
2748
  context_lex_free,
2749
};
2750
2751
static struct isl_tab *context_tab_for_lexmin(struct isl_basic_set *bset)
2752
0
{
2753
0
  struct isl_tab *tab;
2754
0
2755
0
  if (!bset)
2756
0
    return NULL;
2757
0
  tab = tab_for_lexmin(bset_to_bmap(bset), NULL, 1, 0);
2758
0
  if (isl_tab_track_bset(tab, bset) < 0)
2759
0
    goto error;
2760
0
  tab = isl_tab_init_samples(tab);
2761
0
  return tab;
2762
0
error:
2763
0
  isl_tab_free(tab);
2764
0
  return NULL;
2765
0
}
2766
2767
static struct isl_context *isl_context_lex_alloc(struct isl_basic_set *dom)
2768
0
{
2769
0
  struct isl_context_lex *clex;
2770
0
2771
0
  if (!dom)
2772
0
    return NULL;
2773
0
2774
0
  
clex = 0
isl_alloc_type0
(dom->ctx, struct isl_context_lex);
2775
0
  if (!clex)
2776
0
    return NULL;
2777
0
2778
0
  clex->context.op = &isl_context_lex_op;
2779
0
2780
0
  clex->tab = context_tab_for_lexmin(isl_basic_set_copy(dom));
2781
0
  if (restore_lexmin(clex->tab) < 0)
2782
0
    goto error;
2783
0
  clex->tab = check_integer_feasible(clex->tab);
2784
0
  if (!clex->tab)
2785
0
    goto error;
2786
0
2787
0
  return &clex->context;
2788
0
error:
2789
0
  clex->context.op->free(&clex->context);
2790
0
  return NULL;
2791
0
}
2792
2793
/* Representation of the context when using generalized basis reduction.
2794
 *
2795
 * "shifted" contains the offsets of the unit hypercubes that lie inside the
2796
 * context.  Any rational point in "shifted" can therefore be rounded
2797
 * up to an integer point in the context.
2798
 * If the context is constrained by any equality, then "shifted" is not used
2799
 * as it would be empty.
2800
 */
2801
struct isl_context_gbr {
2802
  struct isl_context context;
2803
  struct isl_tab *tab;
2804
  struct isl_tab *shifted;
2805
  struct isl_tab *cone;
2806
};
2807
2808
static struct isl_tab *context_gbr_detect_nonnegative_parameters(
2809
  struct isl_context *context, struct isl_tab *tab)
2810
2.80k
{
2811
2.80k
  struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
2812
2.80k
  if (!tab)
2813
0
    return NULL;
2814
2.80k
  return tab_detect_nonnegative_parameters(tab, cgbr->tab);
2815
2.80k
}
2816
2817
static struct isl_basic_set *context_gbr_peek_basic_set(
2818
  struct isl_context *context)
2819
10.2k
{
2820
10.2k
  struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
2821
10.2k
  if (!cgbr->tab)
2822
0
    return NULL;
2823
10.2k
  return isl_tab_peek_bset(cgbr->tab);
2824
10.2k
}
2825
2826
static struct isl_tab *context_gbr_peek_tab(struct isl_context *context)
2827
13.0k
{
2828
13.0k
  struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
2829
13.0k
  return cgbr->tab;
2830
13.0k
}
2831
2832
/* Initialize the "shifted" tableau of the context, which
2833
 * contains the constraints of the original tableau shifted
2834
 * by the sum of all negative coefficients.  This ensures
2835
 * that any rational point in the shifted tableau can
2836
 * be rounded up to yield an integer point in the original tableau.
2837
 */
2838
static void gbr_init_shifted(struct isl_context_gbr *cgbr)
2839
48
{
2840
48
  int i, j;
2841
48
  struct isl_vec *cst;
2842
48
  struct isl_basic_set *bset = isl_tab_peek_bset(cgbr->tab);
2843
48
  unsigned dim = isl_basic_set_total_dim(bset);
2844
48
2845
48
  cst = isl_vec_alloc(cgbr->tab->mat->ctx, bset->n_ineq);
2846
48
  if (!cst)
2847
0
    return;
2848
48
2849
571
  
for (i = 0; 48
i < bset->n_ineq571
;
++i523
)
{523
2850
523
    isl_int_set(cst->el[i], bset->ineq[i][0]);
2851
5.39k
    for (j = 0; 
j < dim5.39k
;
++j4.87k
)
{4.87k
2852
4.87k
      if (
!4.87k
isl_int_is_neg4.87k
(bset->ineq[i][1 + j]))
2853
4.40k
        continue;
2854
472
      
isl_int_add472
(bset->ineq[i][0], bset->ineq[i][0],472
2855
472
            bset->ineq[i][1 + j]);
2856
472
    }
2857
523
  }
2858
48
2859
48
  cgbr->shifted = isl_tab_from_basic_set(bset, 0);
2860
48
2861
571
  for (i = 0; 
i < bset->n_ineq571
;
++i523
)
2862
523
    isl_int_set(bset->ineq[i][0], cst->el[i]);
2863
48
2864
48
  isl_vec_free(cst);
2865
48
}
2866
2867
/* Check if the shifted tableau is non-empty, and if so
2868
 * use the sample point to construct an integer point
2869
 * of the context tableau.
2870
 */
2871
static struct isl_vec *gbr_get_shifted_sample(struct isl_context_gbr *cgbr)
2872
52
{
2873
52
  struct isl_vec *sample;
2874
52
2875
52
  if (!cgbr->shifted)
2876
48
    gbr_init_shifted(cgbr);
2877
52
  if (!cgbr->shifted)
2878
0
    return NULL;
2879
52
  
if (52
cgbr->shifted->empty52
)
2880
32
    return isl_vec_alloc(cgbr->tab->mat->ctx, 0);
2881
52
2882
20
  sample = isl_tab_get_sample_value(cgbr->shifted);
2883
20
  sample = isl_vec_ceil(sample);
2884
20
2885
20
  return sample;
2886
52
}
2887
2888
static struct isl_basic_set *drop_constant_terms(struct isl_basic_set *bset)
2889
293
{
2890
293
  int i;
2891
293
2892
293
  if (!bset)
2893
0
    return NULL;
2894
293
2895
707
  
for (i = 0; 293
i < bset->n_eq707
;
++i414
)
2896
414
    isl_int_set_si(bset->eq[i][0], 0);
2897
293
2898
2.61k
  for (i = 0; 
i < bset->n_ineq2.61k
;
++i2.32k
)
2899
2.32k
    isl_int_set_si(bset->ineq[i][0], 0);
2900
293
2901
293
  return bset;
2902
293
}
2903
2904
static int use_shifted(struct isl_context_gbr *cgbr)
2905
358
{
2906
358
  if (!cgbr->tab)
2907
0
    return 0;
2908
358
  
return cgbr->tab->bmap->n_eq == 0 && 358
cgbr->tab->bmap->n_div == 0152
;
2909
358
}
2910
2911
static struct isl_vec *gbr_get_sample(struct isl_context_gbr *cgbr)
2912
4.58k
{
2913
4.58k
  struct isl_basic_set *bset;
2914
4.58k
  struct isl_basic_set *cone;
2915
4.58k
2916
4.58k
  if (isl_tab_sample_is_integer(cgbr->tab))
2917
4.23k
    return isl_tab_get_sample_value(cgbr->tab);
2918
4.58k
2919
355
  
if (355
use_shifted(cgbr)355
)
{52
2920
52
    struct isl_vec *sample;
2921
52
2922
52
    sample = gbr_get_shifted_sample(cgbr);
2923
52
    if (
!sample || 52
sample->size > 052
)
2924
20
      return sample;
2925
52
2926
32
    isl_vec_free(sample);
2927
32
  }
2928
355
2929
335
  
if (335
!cgbr->cone335
)
{229
2930
229
    bset = isl_tab_peek_bset(cgbr->tab);
2931
229
    cgbr->cone = isl_tab_from_recession_cone(bset, 0);
2932
229
    if (!cgbr->cone)
2933
0
      return NULL;
2934
229
    
if (229
isl_tab_track_bset(cgbr->cone,229
2935
229
          isl_basic_set_copy(bset)) < 0)
2936
0
      return NULL;
2937
229
  }
2938
335
  
if (335
isl_tab_detect_implicit_equalities(cgbr->cone) < 0335
)
2939
0
    return NULL;
2940
335
2941
335
  
if (335
cgbr->cone->n_dead == cgbr->cone->n_col335
)
{42
2942
42
    struct isl_vec *sample;
2943
42
    struct isl_tab_undo *snap;
2944
42
2945
42
    if (
cgbr->tab->basis42
)
{22
2946
22
      if (
cgbr->tab->basis->n_col != 1 + cgbr->tab->n_var22
)
{6
2947
6
        isl_mat_free(cgbr->tab->basis);
2948
6
        cgbr->tab->basis = NULL;
2949
6
      }
2950
22
      cgbr->tab->n_zero = 0;
2951
22
      cgbr->tab->n_unbounded = 0;
2952
22
    }
2953
42
2954
42
    snap = isl_tab_snap(cgbr->tab);
2955
42
2956
42
    sample = isl_tab_sample(cgbr->tab);
2957
42
2958
42
    if (
!sample || 42
isl_tab_rollback(cgbr->tab, snap) < 042
)
{0
2959
0
      isl_vec_free(sample);
2960
0
      return NULL;
2961
0
    }
2962
42
2963
42
    return sample;
2964
42
  }
2965
335
2966
293
  cone = isl_basic_set_dup(isl_tab_peek_bset(cgbr->cone));
2967
293
  cone = drop_constant_terms(cone);
2968
293
  cone = isl_basic_set_update_from_tab(cone, cgbr->cone);
2969
293
  cone = isl_basic_set_underlying_set(cone);
2970
293
  cone = isl_basic_set_gauss(cone, NULL);
2971
293
2972
293
  bset = isl_basic_set_dup(isl_tab_peek_bset(cgbr->tab));
2973
293
  bset = isl_basic_set_update_from_tab(bset, cgbr->tab);
2974
293
  bset = isl_basic_set_underlying_set(bset);
2975
293
  bset = isl_basic_set_gauss(bset, NULL);
2976
293
2977
293
  return isl_basic_set_sample_with_cone(bset, cone);
2978
335
}
2979
2980
static void check_gbr_integer_feasible(struct isl_context_gbr *cgbr)
2981
15.3k
{
2982
15.3k
  struct isl_vec *sample;
2983
15.3k
2984
15.3k
  if (!cgbr->tab)
2985
0
    return;
2986
15.3k
2987
15.3k
  
if (15.3k
cgbr->tab->empty15.3k
)
2988
10.8k
    return;
2989
15.3k
2990
4.58k
  sample = gbr_get_sample(cgbr);
2991
4.58k
  if (!sample)
2992
0
    goto error;
2993
4.58k
2994
4.58k
  
if (4.58k
sample->size == 04.58k
)
{82
2995
82
    isl_vec_free(sample);
2996
82
    if (isl_tab_mark_empty(cgbr->tab) < 0)
2997
0
      goto error;
2998
82
    return;
2999
82
  }
3000
4.58k
3001
4.50k
  
if (4.50k
isl_tab_add_sample(cgbr->tab, sample) < 04.50k
)
3002
0
    goto error;
3003
4.50k
3004
4.50k
  return;
3005
0
error:
3006
0
  isl_tab_free(cgbr->tab);
3007
0
  cgbr->tab = NULL;
3008
0
}
3009
3010
static struct isl_tab *add_gbr_eq(struct isl_tab *tab, isl_int *eq)
3011
3.54k
{
3012
3.54k
  if (!tab)
3013
0
    return NULL;
3014
3.54k
3015
3.54k
  
if (3.54k
isl_tab_extend_cons(tab, 2) < 03.54k
)
3016
0
    goto error;
3017
3.54k
3018
3.54k
  
if (3.54k
isl_tab_add_eq(tab, eq) < 03.54k
)
3019
0
    goto error;
3020
3.54k
3021
3.54k
  return tab;
3022
0
error:
3023
0
  isl_tab_free(tab);
3024
0
  return NULL;
3025
3.54k
}
3026
3027
/* Add the equality described by "eq" to the context.
3028
 * If "check" is set, then we check if the context is empty after
3029
 * adding the equality.
3030
 * If "update" is set, then we check if the samples are still valid.
3031
 *
3032
 * We do not explicitly add shifted copies of the equality to
3033
 * cgbr->shifted since they would conflict with each other.
3034
 * Instead, we directly mark cgbr->shifted empty.
3035
 */
3036
static void context_gbr_add_eq(struct isl_context *context, isl_int *eq,
3037
    int check, int update)
3038
3.54k
{
3039
3.54k
  struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3040
3.54k
3041
3.54k
  cgbr->tab = add_gbr_eq(cgbr->tab, eq);
3042
3.54k
3043
3.54k
  if (
cgbr->shifted && 3.54k
!cgbr->shifted->empty2
&&
use_shifted(cgbr)1
)
{1
3044
1
    if (isl_tab_mark_empty(cgbr->shifted) < 0)
3045
0
      goto error;
3046
1
  }
3047
3.54k
3048
3.54k
  
if (3.54k
cgbr->cone && 3.54k
cgbr->cone->n_col != cgbr->cone->n_dead105
)
{100
3049
100
    if (isl_tab_extend_cons(cgbr->cone, 2) < 0)
3050
0
      goto error;
3051
100
    
if (100
isl_tab_add_eq(cgbr->cone, eq) < 0100
)
3052
0
      goto error;
3053
100
  }
3054
3.54k
3055
3.54k
  
if (3.54k
check3.54k
)
{3.54k
3056
3.54k
    int v = tab_has_valid_sample(cgbr->tab, eq, 1);
3057
3.54k
    if (v < 0)
3058
0
      goto error;
3059
3.54k
    
if (3.54k
!v3.54k
)
3060
55
      check_gbr_integer_feasible(cgbr);
3061
3.54k
  }
3062
3.54k
  
if (3.54k
update3.54k
)
3063
3.54k
    cgbr->tab = check_samples(cgbr->tab, eq, 1);
3064
3.54k
  return;
3065
0
error:
3066
0
  isl_tab_free(cgbr->tab);
3067
0
  cgbr->tab = NULL;
3068
0
}
3069
3070
static void add_gbr_ineq(struct isl_context_gbr *cgbr, isl_int *ineq)
3071
15.2k
{
3072
15.2k
  if (!cgbr->tab)
3073
0
    return;
3074
15.2k
3075
15.2k
  
if (15.2k
isl_tab_extend_cons(cgbr->tab, 1) < 015.2k
)
3076
0
    goto error;
3077
15.2k
3078
15.2k
  
if (15.2k
isl_tab_add_ineq(cgbr->tab, ineq) < 015.2k
)
3079
0
    goto error;
3080
15.2k
3081
15.2k
  
if (15.2k
cgbr->shifted && 15.2k
!cgbr->shifted->empty15
&&
use_shifted(cgbr)2
)
{2
3082
2
    int i;
3083
2
    unsigned dim;
3084
2
    dim = isl_basic_map_total_dim(cgbr->tab->bmap);
3085
2
3086
2
    if (isl_tab_extend_cons(cgbr->shifted, 1) < 0)
3087
0
      goto error;
3088
2
3089
10
    
for (i = 0; 2
i < dim10
;
++i8
)
{8
3090
8
      if (
!8
isl_int_is_neg8
(ineq[1 + i]))
3091
6
        continue;
3092
2
      
isl_int_add2
(ineq[0], ineq[0], ineq[1 + i]);2
3093
2
    }
3094
2
3095
2
    if (isl_tab_add_ineq(cgbr->shifted, ineq) < 0)
3096
0
      goto error;
3097
2
3098
10
    
for (i = 0; 2
i < dim10
;
++i8
)
{8
3099
8
      if (
!8
isl_int_is_neg8
(ineq[1 + i]))
3100
6
        continue;
3101
2
      
isl_int_sub2
(ineq[0], ineq[0], ineq[1 + i]);2
3102
2
    }
3103
2
  }
3104
15.2k
3105
15.2k
  
if (15.2k
cgbr->cone && 15.2k
cgbr->cone->n_col != cgbr->cone->n_dead2.75k
)
{1.77k
3106
1.77k
    if (isl_tab_extend_cons(cgbr->cone, 1) < 0)
3107
0
      goto error;
3108
1.77k
    
if (1.77k
isl_tab_add_ineq(cgbr->cone, ineq) < 01.77k
)
3109
0
      goto error;
3110
1.77k
  }
3111
15.2k
3112
15.2k
  return;
3113
0
error:
3114
0
  isl_tab_free(cgbr->tab);
3115
0
  cgbr->tab = NULL;
3116
0
}
3117
3118
static void context_gbr_add_ineq(struct isl_context *context, isl_int *ineq,
3119
    int check, int update)
3120
10.2k
{
3121
10.2k
  struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3122
10.2k
3123
10.2k
  add_gbr_ineq(cgbr, ineq);
3124
10.2k
  if (!cgbr->tab)
3125
0
    return;
3126
10.2k
3127
10.2k
  
if (10.2k
check10.2k
)
{7.16k
3128
7.16k
    int v = tab_has_valid_sample(cgbr->tab, ineq, 0);
3129
7.16k
    if (v < 0)
3130
0
      goto error;
3131
7.16k
    
if (7.16k
!v7.16k
)
3132
7.07k
      check_gbr_integer_feasible(cgbr);
3133
7.16k
  }
3134
10.2k
  
if (10.2k
update10.2k
)
3135
2.33k
    cgbr->tab = check_samples(cgbr->tab, ineq, 0);
3136
10.2k
  return;
3137
0
error:
3138
0
  isl_tab_free(cgbr->tab);
3139
0
  cgbr->tab = NULL;
3140
0
}
3141
3142
static isl_stat context_gbr_add_ineq_wrap(void *user, isl_int *ineq)
3143
772
{
3144
772
  struct isl_context *context = (struct isl_context *)user;
3145
772
  context_gbr_add_ineq(context, ineq, 0, 0);
3146
772
  return context->op->is_ok(context) ? 
isl_stat_ok772
:
isl_stat_error0
;
3147
772
}
3148
3149
static enum isl_tab_row_sign context_gbr_ineq_sign(struct isl_context *context,
3150
      isl_int *ineq, int strict)
3151
5.43k
{
3152
5.43k
  struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3153
5.43k
  return tab_ineq_sign(cgbr->tab, ineq, strict);
3154
5.43k
}
3155
3156
/* Check whether "ineq" can be added to the tableau without rendering
3157
 * it infeasible.
3158
 */
3159
static int context_gbr_test_ineq(struct isl_context *context, isl_int *ineq)
3160
5.04k
{
3161
5.04k
  struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3162
5.04k
  struct isl_tab_undo *snap;
3163
5.04k
  struct isl_tab_undo *shifted_snap = NULL;
3164
5.04k
  struct isl_tab_undo *cone_snap = NULL;
3165
5.04k
  int feasible;
3166
5.04k
3167
5.04k
  if (!cgbr->tab)
3168
0
    return -1;
3169
5.04k
3170
5.04k
  
if (5.04k
isl_tab_extend_cons(cgbr->tab, 1) < 05.04k
)
3171
0
    return -1;
3172
5.04k
3173
5.04k
  snap = isl_tab_snap(cgbr->tab);
3174
5.04k
  if (cgbr->shifted)
3175
9
    shifted_snap = isl_tab_snap(cgbr->shifted);
3176
5.04k
  if (cgbr->cone)
3177
1.68k
    cone_snap = isl_tab_snap(cgbr->cone);
3178
5.04k
  add_gbr_ineq(cgbr, ineq);
3179
5.04k
  check_gbr_integer_feasible(cgbr);
3180
5.04k
  if (!cgbr->tab)
3181
0
    return -1;
3182
5.04k
  feasible = !cgbr->tab->empty;
3183
5.04k
  if (isl_tab_rollback(cgbr->tab, snap) < 0)
3184
0
    return -1;
3185
5.04k
  
if (5.04k
shifted_snap5.04k
)
{9
3186
9
    if (isl_tab_rollback(cgbr->shifted, shifted_snap))
3187
0
      return -1;
3188
5.03k
  } else 
if (5.03k
cgbr->shifted5.03k
)
{46
3189
46
    isl_tab_free(cgbr->shifted);
3190
46
    cgbr->shifted = NULL;
3191
46
  }
3192
5.04k
  
if (5.04k
cone_snap5.04k
)
{1.68k
3193
1.68k
    if (isl_tab_rollback(cgbr->cone, cone_snap))
3194
0
      return -1;
3195
3.36k
  } else 
if (3.36k
cgbr->cone3.36k
)
{107
3196
107
    isl_tab_free(cgbr->cone);
3197
107
    cgbr->cone = NULL;
3198
107
  }
3199
5.04k
3200
5.04k
  return feasible;
3201
5.04k
}
3202
3203
/* Return the column of the last of the variables associated to
3204
 * a column that has a non-zero coefficient.
3205
 * This function is called in a context where only coefficients
3206
 * of parameters or divs can be non-zero.
3207
 */
3208
static int last_non_zero_var_col(struct isl_tab *tab, isl_int *p)
3209
364
{
3210
364
  int i;
3211
364
  int col;
3212
364
3213
364
  if (tab->n_var == 0)
3214
0
    return -1;
3215
364
3216
792
  
for (i = tab->n_var - 1; 364
i >= 0792
;
--i428
)
{792
3217
792
    if (
i >= tab->n_param && 792
i < tab->n_var - tab->n_div779
)
3218
24
      continue;
3219
768
    
if (768
tab->var[i].is_row768
)
3220
115
      continue;
3221
653
    col = tab->var[i].index;
3222
653
    if (
!653
isl_int_is_zero653
(p[col]))
3223
364
      return col;
3224
653
  }
3225
364
3226
0
  return -1;
3227
364
}
3228
3229
/* Look through all the recently added equalities in the context
3230
 * to see if we can propagate any of them to the main tableau.
3231
 *
3232
 * The newly added equalities in the context are encoded as pairs
3233
 * of inequalities starting at inequality "first".
3234
 *
3235
 * We tentatively add each of these equalities to the main tableau
3236
 * and if this happens to result in a row with a final coefficient
3237
 * that is one or negative one, we use it to kill a column
3238
 * in the main tableau.  Otherwise, we discard the tentatively
3239
 * added row.
3240
 * This tentative addition of equality constraints turns
3241
 * on the undo facility of the tableau.  Turn it off again
3242
 * at the end, assuming it was turned off to begin with.
3243
 *
3244
 * Return 0 on success and -1 on failure.
3245
 */
3246
static int propagate_equalities(struct isl_context_gbr *cgbr,
3247
  struct isl_tab *tab, unsigned first)
3248
184
{
3249
184
  int i;
3250
184
  struct isl_vec *eq = NULL;
3251
184
  isl_bool needs_undo;
3252
184
3253
184
  needs_undo = isl_tab_need_undo(tab);
3254
184
  if (needs_undo < 0)
3255
0
    goto error;
3256
184
  eq = isl_vec_alloc(tab->mat->ctx, 1 + tab->n_var);
3257
184
  if (!eq)
3258
0
    goto error;
3259
184
3260
184
  
if (184
isl_tab_extend_cons(tab, (cgbr->tab->bmap->n_ineq - first)/2) < 0184
)
3261
0
    goto error;
3262
184
3263
184
  isl_seq_clr(eq->el + 1 + tab->n_param,
3264
184
        tab->n_var - tab->n_param - tab->n_div);
3265
548
  for (i = first; 
i < cgbr->tab->bmap->n_ineq548
;
i += 2364
)
{364
3266
364
    int j;
3267
364
    int r;
3268
364
    struct isl_tab_undo *snap;
3269
364
    snap = isl_tab_snap(tab);
3270
364
3271
364
    isl_seq_cpy(eq->el, cgbr->tab->bmap->ineq[i], 1 + tab->n_param);
3272
364
    isl_seq_cpy(eq->el + 1 + tab->n_var - tab->n_div,
3273
364
          cgbr->tab->bmap->ineq[i] + 1 + tab->n_param,
3274
364
          tab->n_div);
3275
364
3276
364
    r = isl_tab_add_row(tab, eq->el);
3277
364
    if (r < 0)
3278
0
      goto error;
3279
364
    r = tab->con[r].index;
3280
364
    j = last_non_zero_var_col(tab, tab->mat->row[r] + 2 + tab->M);
3281
364
    if (
j < 0 || 364
j < tab->n_dead364
||
3282
364
        
!364
isl_int_is_one364
(tab->mat->row[r][0]) ||
3283
364
        
(!364
isl_int_is_one364
(tab->mat->row[r][2 + tab->M + j]) &&
3284
314
         
!314
isl_int_is_negone314
(tab->mat->row[r][2 + tab->M + j])))
{260
3285
260
      if (isl_tab_rollback(tab, snap) < 0)
3286
0
        goto error;
3287
260
      continue;
3288
260
    }
3289
104
    
if (104
isl_tab_pivot(tab, r, j) < 0104
)
3290
0
      goto error;
3291
104
    
if (104
isl_tab_kill_col(tab, j) < 0104
)
3292
0
      goto error;
3293
104
3294
104
    
if (104
restore_lexmin(tab) < 0104
)
3295
0
      goto error;
3296
104
  }
3297
184
3298
184
  
if (184
!needs_undo184
)
3299
184
    isl_tab_clear_undo(tab);
3300
184
  isl_vec_free(eq);
3301
184
3302
184
  return 0;
3303
0
error:
3304
0
  isl_vec_free(eq);
3305
0
  isl_tab_free(cgbr->tab);
3306
0
  cgbr->tab = NULL;
3307
0
  return -1;
3308
184
}
3309
3310
static int context_gbr_detect_equalities(struct isl_context *context,
3311
  struct isl_tab *tab)
3312
363
{
3313
363
  struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3314
363
  unsigned n_ineq;
3315
363
3316
363
  if (
!cgbr->cone363
)
{96
3317
96
    struct isl_basic_set *bset = isl_tab_peek_bset(cgbr->tab);
3318
96
    cgbr->cone = isl_tab_from_recession_cone(bset, 0);
3319
96
    if (!cgbr->cone)
3320
0
      goto error;
3321
96
    
if (96
isl_tab_track_bset(cgbr->cone,96
3322
96
          isl_basic_set_copy(bset)) < 0)
3323
0
      goto error;
3324
96
  }
3325
363
  
if (363
isl_tab_detect_implicit_equalities(cgbr->cone) < 0363
)
3326
0
    goto error;
3327
363
3328
363
  n_ineq = cgbr->tab->bmap->n_ineq;
3329
363
  cgbr->tab = isl_tab_detect_equalities(cgbr->tab, cgbr->cone);
3330
363
  if (!cgbr->tab)
3331
0
    return -1;
3332
363
  
if (363
cgbr->tab->bmap->n_ineq > n_ineq &&363
3333
184
      propagate_equalities(cgbr, tab, n_ineq) < 0)
3334
0
    return -1;
3335
363
3336
363
  return 0;
3337
0
error:
3338
0
  isl_tab_free(cgbr->tab);
3339
0
  cgbr->tab = NULL;
3340
0
  return -1;
3341
363
}
3342
3343
static int context_gbr_get_div(struct isl_context *context, struct isl_tab *tab,
3344
    struct isl_vec *div)
3345
401
{
3346
401
  return get_div(tab, context, div);
3347
401
}
3348
3349
static isl_bool context_gbr_insert_div(struct isl_context *context, int pos,
3350
  __isl_keep isl_vec *div)
3351
386
{
3352
386
  struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3353
386
  if (
cgbr->cone386
)
{270
3354
270
    int r, n_div, o_div;
3355
270
3356
270
    n_div = isl_basic_map_dim(cgbr->cone->bmap, isl_dim_div);
3357
270
    o_div = cgbr->cone->n_var - n_div;
3358
270
3359
270
    if (isl_tab_extend_cons(cgbr->cone, 3) < 0)
3360
0
      return isl_bool_error;
3361
270
    
if (270
isl_tab_extend_vars(cgbr->cone, 1) < 0270
)
3362
0
      return isl_bool_error;
3363
270
    
if (270
(r = isl_tab_insert_var(cgbr->cone, pos)) <0270
)
3364
0
      return isl_bool_error;
3365
270
3366
270
    cgbr->cone->bmap = isl_basic_map_insert_div(cgbr->cone->bmap,
3367
270
                r - o_div, div);
3368
270
    if (!cgbr->cone->bmap)
3369
0
      return isl_bool_error;
3370
270
    
if (270
isl_tab_push_var(cgbr->cone, isl_tab_undo_bmap_div,270
3371
270
            &cgbr->cone->var[r]) < 0)
3372
0
      return isl_bool_error;
3373
270
  }
3374
386
  return context_tab_insert_div(cgbr->tab, pos, div,
3375
386
          context_gbr_add_ineq_wrap, context);
3376
386
}
3377
3378
static int context_gbr_best_split(struct isl_context *context,
3379
    struct isl_tab *tab)
3380
235
{
3381
235
  struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3382
235
  struct isl_tab_undo *snap;
3383
235
  int r;
3384
235
3385
235
  snap = isl_tab_snap(cgbr->tab);
3386
235
  r = best_split(tab, cgbr->tab);
3387
235
3388
235
  if (
r >= 0 && 235
isl_tab_rollback(cgbr->tab, snap) < 0235
)
3389
0
    return -1;
3390
235
3391
235
  return r;
3392
235
}
3393
3394
static int context_gbr_is_empty(struct isl_context *context)
3395
18.0k
{
3396
18.0k
  struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3397
18.0k
  if (!cgbr->tab)
3398
0
    return -1;
3399
18.0k
  return cgbr->tab->empty;
3400
18.0k
}
3401
3402
struct isl_gbr_tab_undo {
3403
  struct isl_tab_undo *tab_snap;
3404
  struct isl_tab_undo *shifted_snap;
3405
  struct isl_tab_undo *cone_snap;
3406
};
3407
3408
static void *context_gbr_save(struct isl_context *context)
3409
11.0k
{
3410
11.0k
  struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3411
11.0k
  struct isl_gbr_tab_undo *snap;
3412
11.0k
3413
11.0k
  if (!cgbr->tab)
3414
0
    return NULL;
3415
11.0k
3416
11.0k
  
snap = 11.0k
isl_alloc_type11.0k
(cgbr->tab->mat->ctx, struct isl_gbr_tab_undo);
3417
11.0k
  if (!snap)
3418
0
    return NULL;
3419
11.0k
3420
11.0k
  snap->tab_snap = isl_tab_snap(cgbr->tab);
3421
11.0k
  if (isl_tab_save_samples(cgbr->tab) < 0)
3422
0
    goto error;
3423
11.0k
3424
11.0k
  
if (11.0k
cgbr->shifted11.0k
)
3425
6
    snap->shifted_snap = isl_tab_snap(cgbr->shifted);
3426
11.0k
  else
3427
11.0k
    snap->shifted_snap = NULL;
3428
11.0k
3429
11.0k
  if (cgbr->cone)
3430
405
    snap->cone_snap = isl_tab_snap(cgbr->cone);
3431
11.0k
  else
3432
10.6k
    snap->cone_snap = NULL;
3433
11.0k
3434
11.0k
  return snap;
3435
0
error:
3436
0
  free(snap);
3437
0
  return NULL;
3438
11.0k
}
3439
3440
static void context_gbr_restore(struct isl_context *context, void *save)
3441
8.89k
{
3442
8.89k
  struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3443
8.89k
  struct isl_gbr_tab_undo *snap = (struct isl_gbr_tab_undo *)save;
3444
8.89k
  if (!snap)
3445
0
    goto error;
3446
8.89k
  
if (8.89k
isl_tab_rollback(cgbr->tab, snap->tab_snap) < 08.89k
)
3447
0
    goto error;
3448
8.89k
3449
8.89k
  
if (8.89k
snap->shifted_snap8.89k
)
{4
3450
4
    if (isl_tab_rollback(cgbr->shifted, snap->shifted_snap) < 0)
3451
0
      goto error;
3452
8.89k
  } else 
if (8.89k
cgbr->shifted8.89k
)
{0
3453
0
    isl_tab_free(cgbr->shifted);
3454
0
    cgbr->shifted = NULL;
3455
0
  }
3456
8.89k
3457
8.89k
  
if (8.89k
snap->cone_snap8.89k
)
{399
3458
399
    if (isl_tab_rollback(cgbr->cone, snap->cone_snap) < 0)
3459
0
      goto error;
3460
8.49k
  } else 
if (8.49k
cgbr->cone8.49k
)
{105
3461
105
    isl_tab_free(cgbr->cone);
3462
105
    cgbr->cone = NULL;
3463
105
  }
3464
8.89k
3465
8.89k
  free(snap);
3466
8.89k
3467
8.89k
  return;
3468
0
error:
3469
0
  free(snap);
3470
0
  isl_tab_free(cgbr->tab);
3471
0
  cgbr->tab = NULL;
3472
0
}
3473
3474
static void context_gbr_discard(void *save)
3475
2.18k
{
3476
2.18k
  struct isl_gbr_tab_undo *snap = (struct isl_gbr_tab_undo *)save;
3477
2.18k
  free(snap);
3478
2.18k
}
3479
3480
static int context_gbr_is_ok(struct isl_context *context)
3481
1.19k
{
3482
1.19k
  struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3483
1.19k
  return !!cgbr->tab;
3484
1.19k
}
3485
3486
static void context_gbr_invalidate(struct isl_context *context)
3487
0
{
3488
0
  struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3489
0
  isl_tab_free(cgbr->tab);
3490
0
  cgbr->tab = NULL;
3491
0
}
3492
3493
static __isl_null struct isl_context *context_gbr_free(
3494
  struct isl_context *context)
3495
3.22k
{
3496
3.22k
  struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3497
3.22k
  isl_tab_free(cgbr->tab);
3498
3.22k
  isl_tab_free(cgbr->shifted);
3499
3.22k
  isl_tab_free(cgbr->cone);
3500
3.22k
  free(cgbr);
3501
3.22k
3502
3.22k
  return NULL;
3503
3.22k
}
3504
3505
struct isl_context_op isl_context_gbr_op = {
3506
  context_gbr_detect_nonnegative_parameters,
3507
  context_gbr_peek_basic_set,
3508
  context_gbr_peek_tab,
3509
  context_gbr_add_eq,
3510
  context_gbr_add_ineq,
3511
  context_gbr_ineq_sign,
3512
  context_gbr_test_ineq,
3513
  context_gbr_get_div,
3514
  context_gbr_insert_div,
3515
  context_gbr_detect_equalities,
3516
  context_gbr_best_split,
3517
  context_gbr_is_empty,
3518
  context_gbr_is_ok,
3519
  context_gbr_save,
3520
  context_gbr_restore,
3521
  context_gbr_discard,
3522
  context_gbr_invalidate,
3523
  context_gbr_free,
3524
};
3525
3526
static struct isl_context *isl_context_gbr_alloc(__isl_keep isl_basic_set *dom)
3527
3.22k
{
3528
3.22k
  struct isl_context_gbr *cgbr;
3529
3.22k
3530
3.22k
  if (!dom)
3531
0
    return NULL;
3532
3.22k
3533
3.22k
  
cgbr = 3.22k
isl_calloc_type3.22k
(dom->ctx, struct isl_context_gbr);
3534
3.22k
  if (!cgbr)
3535
0
    return NULL;
3536
3.22k
3537
3.22k
  cgbr->context.op = &isl_context_gbr_op;
3538
3.22k
3539
3.22k
  cgbr->shifted = NULL;
3540
3.22k
  cgbr->cone = NULL;
3541
3.22k
  cgbr->tab = isl_tab_from_basic_set(dom, 1);
3542
3.22k
  cgbr->tab = isl_tab_init_samples(cgbr->tab);
3543
3.22k
  if (!cgbr->tab)
3544
0
    goto error;
3545
3.22k
  check_gbr_integer_feasible(cgbr);
3546
3.22k
3547
3.22k
  return &cgbr->context;
3548
0
error:
3549
0
  cgbr->context.op->free(&cgbr->context);
3550
0
  return NULL;
3551
3.22k
}
3552
3553
/* Allocate a context corresponding to "dom".
3554
 * The representation specific fields are initialized by
3555
 * isl_context_lex_alloc or isl_context_gbr_alloc.
3556
 * The shared "n_unknown" field is initialized to the number
3557
 * of final unknown integer divisions in "dom".
3558
 */
3559
static struct isl_context *isl_context_alloc(__isl_keep isl_basic_set *dom)
3560
3.22k
{
3561
3.22k
  struct isl_context *context;
3562
3.22k
  int first;
3563
3.22k
3564
3.22k
  if (!dom)
3565
0
    return NULL;
3566
3.22k
3567
3.22k
  
if (3.22k
dom->ctx->opt->context == 3.22k
ISL_CONTEXT_LEXMIN3.22k
)
3568
0
    context = isl_context_lex_alloc(dom);
3569
3.22k
  else
3570
3.22k
    context = isl_context_gbr_alloc(dom);
3571
3.22k
3572
3.22k
  if (!context)
3573
0
    return NULL;
3574
3.22k
3575
3.22k
  first = isl_basic_set_first_unknown_div(dom);
3576
3.22k
  if (first < 0)
3577
0
    return context->op->free(context);
3578
3.22k
  context->n_unknown = isl_basic_set_dim(dom, isl_dim_div) - first;
3579
3.22k
3580
3.22k
  return context;
3581
3.22k
}
3582
3583
/* Initialize some common fields of "sol", which keeps track
3584
 * of the solution of an optimization problem on "bmap" over
3585
 * the domain "dom".
3586
 * If "max" is set, then a maximization problem is being solved, rather than
3587
 * a minimization problem, which means that the variables in the
3588
 * tableau have value "M - x" rather than "M + x".
3589
 */
3590
static isl_stat sol_init(struct isl_sol *sol, __isl_keep isl_basic_map *bmap,
3591
  __isl_keep isl_basic_set *dom, int max)
3592
3.22k
{
3593
3.22k
  sol->rational = ISL_F_ISSET(bmap, ISL_BASIC_MAP_RATIONAL);
3594
3.22k
  sol->dec_level.callback.run = &sol_dec_level_wrap;
3595
3.22k
  sol->dec_level.sol = sol;
3596
3.22k
  sol->max = max;
3597
3.22k
  sol->n_out = isl_basic_map_dim(bmap, isl_dim_out);
3598
3.22k
  sol->space = isl_basic_map_get_space(bmap);
3599
3.22k
3600
3.22k
  sol->context = isl_context_alloc(dom);
3601
3.22k
  if (
!sol->space || 3.22k
!sol->context3.22k
)
3602
0
    return isl_stat_error;
3603
3.22k
3604
3.22k
  return isl_stat_ok;
3605
3.22k
}
3606
3607
/* Construct an isl_sol_map structure for accumulating the solution.
3608
 * If track_empty is set, then we also keep track of the parts
3609
 * of the context where there is no solution.
3610
 * If max is set, then we are solving a maximization, rather than
3611
 * a minimization problem, which means that the variables in the
3612
 * tableau have value "M - x" rather than "M + x".
3613
 */
3614
static struct isl_sol *sol_map_init(__isl_keep isl_basic_map *bmap,
3615
  __isl_take isl_basic_set *dom, int track_empty, int max)
3616
1.27k
{
3617
1.27k
  struct isl_sol_map *sol_map = NULL;
3618
1.27k
  isl_space *space;
3619
1.27k
3620
1.27k
  if (!bmap)
3621
0
    goto error;
3622
1.27k
3623
1.27k
  
sol_map = 1.27k
isl_calloc_type1.27k
(bmap->ctx, struct isl_sol_map);
3624
1.27k
  if (!sol_map)
3625
0
    goto error;
3626
1.27k
3627
1.27k
  sol_map->sol.free = &sol_map_free;
3628
1.27k
  if (sol_init(&sol_map->sol, bmap, dom, max) < 0)
3629
0
    goto error;
3630
1.27k
  sol_map->sol.add = &sol_map_add_wrap;
3631
941
  sol_map->sol.add_empty = track_empty ? &sol_map_add_empty_wrap : NULL;
3632
1.27k
  space = isl_space_copy(sol_map->sol.space);
3633
1.27k
  sol_map->map = isl_map_alloc_space(space, 1, ISL_MAP_DISJOINT);
3634
1.27k
  if (!sol_map->map)
3635
0
    goto error;
3636
1.27k
3637
1.27k
  
if (1.27k
track_empty1.27k
)
{941
3638
941
    sol_map->empty = isl_set_alloc_space(isl_basic_set_get_space(dom),
3639
941
              1, ISL_SET_DISJOINT);
3640
941
    if (!sol_map->empty)
3641
0
      goto error;
3642
941
  }
3643
1.27k
3644
1.27k
  isl_basic_set_free(dom);
3645
1.27k
  return &sol_map->sol;
3646
0
error:
3647
0
  isl_basic_set_free(dom);
3648
0
  sol_free(&sol_map->sol);
3649
0
  return NULL;
3650
1.27k
}
3651
3652
/* Check whether all coefficients of (non-parameter) variables
3653
 * are non-positive, meaning that no pivots can be performed on the row.
3654
 */
3655
static int is_critical(struct isl_tab *tab, int row)
3656
5.43k
{
3657
5.43k
  int j;
3658
5.43k
  unsigned off = 2 + tab->M;
3659
5.43k
3660
22.9k
  for (j = tab->n_dead; 
j < tab->n_col22.9k
;
++j17.5k
)
{19.0k
3661
19.0k
    if (tab->col_var[j] >= 0 &&
3662
11.8k
        (tab->col_var[j] < tab->n_param  ||
3663
1.16k
        tab->col_var[j] >= tab->n_var - tab->n_div))
3664
11.8k
      continue;
3665
19.0k
3666
7.17k
    
if (7.17k
isl_int_is_pos7.17k
(tab->mat->row[row][off + j]))
3667
1.48k
      return 0;
3668
7.17k
  }
3669
5.43k
3670
3.94k
  return 1;
3671
5.43k
}
3672
3673
/* Check whether the inequality represented by vec is strict over the integers,
3674
 * i.e., there are no integer values satisfying the constraint with
3675
 * equality.  This happens if the gcd of the coefficients is not a divisor
3676
 * of the constant term.  If so, scale the constraint down by the gcd
3677
 * of the coefficients.
3678
 */
3679
static int is_strict(struct isl_vec *vec)
3680
6.58k
{
3681
6.58k
  isl_int gcd;
3682
6.58k
  int strict = 0;
3683
6.58k
3684
6.58k
  isl_int_init(gcd);
3685
6.58k
  isl_seq_gcd(vec->el + 1, vec->size - 1, &gcd);
3686
6.58k
  if (
!6.58k
isl_int_is_one6.58k
(gcd))
{948
3687
948
    strict = !isl_int_is_divisible_by(vec->el[0], gcd);
3688
948
    isl_int_fdiv_q(vec->el[0], vec->el[0], gcd);
3689
948
    isl_seq_scale_down(vec->el + 1, vec->el + 1, gcd, vec->size-1);
3690
948
  }
3691
6.58k
  isl_int_clear(gcd);
3692
6.58k
3693
6.58k
  return strict;
3694
6.58k
}
3695
3696
/* Determine the sign of the given row of the main tableau.
3697
 * The result is one of
3698
 *  isl_tab_row_pos: always non-negative; no pivot needed
3699
 *  isl_tab_row_neg: always non-positive; pivot
3700
 *  isl_tab_row_any: can be both positive and negative; split
3701
 *
3702
 * We first handle some simple cases
3703
 *  - the row sign may be known already
3704
 *  - the row may be obviously non-negative
3705
 *  - the parametric constant may be equal to that of another row
3706
 *    for which we know the sign.  This sign will be either "pos" or
3707
 *    "any".  If it had been "neg" then we would have pivoted before.
3708
 *
3709
 * If none of these cases hold, we check the value of the row for each
3710
 * of the currently active samples.  Based on the signs of these values
3711
 * we make an initial determination of the sign of the row.
3712
 *
3713
 *  all zero      ->  unk(nown)
3714
 *  all non-negative    ->  pos
3715
 *  all non-positive    ->  neg
3716
 *  both negative and positive  ->  all
3717
 *
3718
 * If we end up with "all", we are done.
3719
 * Otherwise, we perform a check for positive and/or negative
3720
 * values as follows.
3721
 *
3722
 *  samples        neg         unk         pos
3723
 *  <0 ?          Y        N      Y        N
3724
 *              pos    any      pos
3725
 *  >0 ?       Y      N  Y     N
3726
 *        any    neg  any   neg
3727
 *
3728
 * There is no special sign for "zero", because we can usually treat zero
3729
 * as either non-negative or non-positive, whatever works out best.
3730
 * However, if the row is "critical", meaning that pivoting is impossible
3731
 * then we don't want to limp zero with the non-positive case, because
3732
 * then we we would lose the solution for those values of the parameters
3733
 * where the value of the row is zero.  Instead, we treat 0 as non-negative
3734
 * ensuring a split if the row can attain both zero and negative values.
3735
 * The same happens when the original constraint was one that could not
3736
 * be satisfied with equality by any integer values of the parameters.
3737
 * In this case, we normalize the constraint, but then a value of zero
3738
 * for the normalized constraint is actually a positive value for the
3739
 * original constraint, so again we need to treat zero as non-negative.
3740
 * In both these cases, we have the following decision tree instead:
3741
 *
3742
 *  all non-negative    ->  pos
3743
 *  all negative      ->  neg
3744
 *  both negative and non-negative  ->  all
3745
 *
3746
 *  samples        neg                     pos
3747
 *  <0 ?                        Y        N
3748
 *                     any      pos
3749
 *  >=0 ?      Y      N
3750
 *        any    neg
3751
 */
3752
static enum isl_tab_row_sign row_sign(struct isl_tab *tab,
3753
  struct isl_sol *sol, int row)
3754
33.6k
{
3755
33.6k
  struct isl_vec *ineq = NULL;
3756
33.6k
  enum isl_tab_row_sign res = isl_tab_row_unknown;
3757
33.6k
  int critical;
3758
33.6k
  int strict;
3759
33.6k
  int row2;
3760
33.6k
3761
33.6k
  if (tab->row_sign[row] != isl_tab_row_unknown)
3762
18.8k
    return tab->row_sign[row];
3763
14.7k
  
if (14.7k
is_obviously_nonneg(tab, row)14.7k
)
3764
9.26k
    return isl_tab_row_pos;
3765
63.2k
  
for (row2 = tab->n_redundant; 5.47k
row2 < tab->n_row63.2k
;
++row257.8k
)
{57.8k
3766
57.8k
    if (tab->row_sign[row2] == isl_tab_row_unknown)
3767
19.0k
      continue;
3768
38.8k
    
if (38.8k
identical_parameter_line(tab, row, row2)38.8k
)
3769
47
      return tab->row_sign[row2];
3770
38.8k
  }
3771
5.47k
3772
5.43k
  critical = is_critical(tab, row);
3773
5.43k
3774
5.43k
  ineq = get_row_parameter_ineq(tab, row);
3775
5.43k
  if (!ineq)
3776
0
    goto error;
3777
5.43k
3778
5.43k
  strict = is_strict(ineq);
3779
5.43k
3780
5.43k
  res = sol->context->op->ineq_sign(sol->context, ineq->el,
3781
1.48k
            critical || strict);
3782
5.43k
3783
5.43k
  if (
res == isl_tab_row_unknown || 5.43k
res == isl_tab_row_pos5.40k
)
{4.11k
3784
4.11k
    /* test for negative values */
3785
4.11k
    int feasible;
3786
4.11k
    isl_seq_neg(ineq->el, ineq->el, ineq->size);
3787
4.11k
    isl_int_sub_ui(ineq->el[0], ineq->el[0], 1);
3788
4.11k
3789
4.11k
    feasible = sol->context->op->test_ineq(sol->context, ineq->el);
3790
4.11k
    if (feasible < 0)
3791
0
      goto error;
3792
4.11k
    
if (4.11k
!feasible4.11k
)
3793
3.83k
      res = isl_tab_row_pos;
3794
4.11k
    else
3795
282
      
res = (res == isl_tab_row_unknown) ? 282
isl_tab_row_neg4
3796
278
                 : isl_tab_row_any;
3797
4.11k
    if (
res == isl_tab_row_neg4.11k
)
{4
3798
4
      isl_seq_neg(ineq->el, ineq->el, ineq->size);
3799
4
      isl_int_sub_ui(ineq->el[0], ineq->el[0], 1);
3800
4
    }
3801
4.11k
  }
3802
5.43k
3803
5.43k
  
if (5.43k
res == isl_tab_row_neg5.43k
)
{935
3804
935
    /* test for positive values */
3805
935
    int feasible;
3806
935
    if (
!critical && 935
!strict146
)
3807
142
      isl_int_sub_ui(ineq->el[0], ineq->el[0], 1);
3808
935
3809
935
    feasible = sol->context->op->test_ineq(sol->context, ineq->el);
3810
935
    if (feasible < 0)
3811
0
      goto error;
3812
935
    
if (935
feasible935
)
3813
836
      res = isl_tab_row_any;
3814
935
  }
3815
5.43k
3816
5.43k
  isl_vec_free(ineq);
3817
5.43k
  return res;
3818
0
error:
3819
0
  isl_vec_free(ineq);
3820
0
  return isl_tab_row_unknown;
3821
5.43k
}
3822
3823
static void find_solutions(struct isl_sol *sol, struct isl_tab *tab);
3824
3825
/* Find solutions for values of the parameters that satisfy the given
3826
 * inequality.
3827
 *
3828
 * We currently take a snapshot of the context tableau that is reset
3829
 * when we return from this function, while we make a copy of the main
3830
 * tableau, leaving the original main tableau untouched.
3831
 * These are fairly arbitrary choices.  Making a copy also of the context
3832
 * tableau would obviate the need to undo any changes made to it later,
3833
 * while taking a snapshot of the main tableau could reduce memory usage.
3834
 * If we were to switch to taking a snapshot of the main tableau,
3835
 * we would have to keep in mind that we need to save the row signs
3836
 * and that we need to do this before saving the current basis
3837
 * such that the basis has been restore before we restore the row signs.
3838
 */
3839
static void find_in_pos(struct isl_sol *sol, struct isl_tab *tab, isl_int *ineq)
3840
1.14k
{
3841
1.14k
  void *saved;
3842
1.14k
3843
1.14k
  if (!sol->context)
3844
0
    goto error;
3845
1.14k
  saved = sol->context->op->save(sol->context);
3846
1.14k
3847
1.14k
  tab = isl_tab_dup(tab);
3848
1.14k
  if (!tab)
3849
0
    goto error;
3850
1.14k
3851
1.14k
  sol->context->op->add_ineq(sol->context, ineq, 0, 1);
3852
1.14k
3853
1.14k
  find_solutions(sol, tab);
3854
1.14k
3855
1.14k
  if (!sol->error)
3856
1.14k
    sol->context->op->restore(sol->context, saved);
3857
1.14k
  else
3858
0
    sol->context->op->discard(saved);
3859
1.14k
  return;
3860
0
error:
3861
0
  sol->error = 1;
3862
0
}
3863
3864
/* Record the absence of solutions for those values of the parameters
3865
 * that do not satisfy the given inequality with equality.
3866
 */
3867
static void no_sol_in_strict(struct isl_sol *sol,
3868
  struct isl_tab *tab, struct isl_vec *ineq)
3869
7.12k
{
3870
7.12k
  int empty;
3871
7.12k
  void *saved;
3872
7.12k
3873
7.12k
  if (
!sol->context || 7.12k
sol->error7.12k
)
3874
0
    goto error;
3875
7.12k
  saved = sol->context->op->save(sol->context);
3876
7.12k
3877
7.12k
  isl_int_sub_ui(ineq->el[0], ineq->el[0], 1);
3878
7.12k
3879
7.12k
  sol->context->op->add_ineq(sol->context, ineq->el, 1, 0);
3880
7.12k
  if (!sol->context)
3881
0
    goto error;
3882
7.12k
3883
7.12k
  empty = tab->empty;
3884
7.12k
  tab->empty = 1;
3885
7.12k
  sol_add(sol, tab);
3886
7.12k
  tab->empty = empty;
3887
7.12k
3888
7.12k
  isl_int_add_ui(ineq->el[0], ineq->el[0], 1);
3889
7.12k
3890
7.12k
  sol->context->op->restore(sol->context, saved);
3891
7.12k
  return;
3892
0
error:
3893
0
  sol->error = 1;
3894
0
}
3895
3896
/* Reset all row variables that are marked to have a sign that may
3897
 * be both positive and negative to have an unknown sign.
3898
 */
3899
static void reset_any_to_unknown(struct isl_tab *tab)
3900
1.14k
{
3901
1.14k
  int row;
3902
1.14k
3903
9.93k
  for (row = tab->n_redundant; 
row < tab->n_row9.93k
;
++row8.78k
)
{8.78k
3904
8.78k
    if (!isl_tab_var_from_row(tab, row)->is_nonneg)
3905
108
      continue;
3906
8.68k
    
if (8.68k
tab->row_sign[row] == isl_tab_row_any8.68k
)
3907
1.48k
      tab->row_sign[row] = isl_tab_row_unknown;
3908
8.68k
  }
3909
1.14k
}
3910
3911
/* Compute the lexicographic minimum of the set represented by the main
3912
 * tableau "tab" within the context "sol->context_tab".
3913
 * On entry the sample value of the main tableau is lexicographically
3914
 * less than or equal to this lexicographic minimum.
3915
 * Pivots are performed until a feasible point is found, which is then
3916
 * necessarily equal to the minimum, or until the tableau is found to
3917
 * be infeasible.  Some pivots may need to be performed for only some
3918
 * feasible values of the context tableau.  If so, the context tableau
3919
 * is split into a part where the pivot is needed and a part where it is not.
3920
 *
3921
 * Whenever we enter the main loop, the main tableau is such that no
3922
 * "obvious" pivots need to be performed on it, where "obvious" means
3923
 * that the given row can be seen to be negative without looking at
3924
 * the context tableau.  In particular, for non-parametric problems,
3925
 * no pivots need to be performed on the main tableau.
3926
 * The caller of find_solutions is responsible for making this property
3927
 * hold prior to the first iteration of the loop, while restore_lexmin
3928
 * is called before every other iteration.
3929
 *
3930
 * Inside the main loop, we first examine the signs of the rows of
3931
 * the main tableau within the context of the context tableau.
3932
 * If we find a row that is always non-positive for all values of
3933
 * the parameters satisfying the context tableau and negative for at
3934
 * least one value of the parameters, we perform the appropriate pivot
3935
 * and start over.  An exception is the case where no pivot can be
3936
 * performed on the row.  In this case, we require that the sign of
3937
 * the row is negative for all values of the parameters (rather than just
3938
 * non-positive).  This special case is handled inside row_sign, which
3939
 * will say that the row can have any sign if it determines that it can
3940
 * attain both negative and zero values.
3941
 *
3942
 * If we can't find a row that always requires a pivot, but we can find
3943
 * one or more rows that require a pivot for some values of the parameters
3944
 * (i.e., the row can attain both positive and negative signs), then we split
3945
 * the context tableau into two parts, one where we force the sign to be
3946
 * non-negative and one where we force is to be negative.
3947
 * The non-negative part is handled by a recursive call (through find_in_pos).
3948
 * Upon returning from this call, we continue with the negative part and
3949
 * perform the required pivot.
3950
 *
3951
 * If no such rows can be found, all rows are non-negative and we have
3952
 * found a (rational) feasible point.  If we only wanted a rational point
3953
 * then we are done.
3954
 * Otherwise, we check if all values of the sample point of the tableau
3955
 * are integral for the variables.  If so, we have found the minimal
3956
 * integral point and we are done.
3957
 * If the sample point is not integral, then we need to make a distinction
3958
 * based on whether the constant term is non-integral or the coefficients
3959
 * of the parameters.  Furthermore, in order to decide how to handle
3960
 * the non-integrality, we also need to know whether the coefficients
3961
 * of the other columns in the tableau are integral.  This leads
3962
 * to the following table.  The first two rows do not correspond
3963
 * to a non-integral sample point and are only mentioned for completeness.
3964
 *
3965
 *  constant  parameters  other
3966
 *
3967
 *  int   int   int |
3968
 *  int   int   rat | -> no problem
3969
 *
3970
 *  rat   int   int   -> fail
3971
 *
3972
 *  rat   int   rat   -> cut
3973
 *
3974
 *  int   rat   rat |
3975
 *  rat   rat   rat | -> parametric cut
3976
 *
3977
 *  int   rat   int |
3978
 *  rat   rat   int | -> split context
3979
 *
3980
 * If the parametric constant is completely integral, then there is nothing
3981
 * to be done.  If the constant term is non-integral, but all the other
3982
 * coefficient are integral, then there is nothing that can be done
3983
 * and the tableau has no integral solution.
3984
 * If, on the other hand, one or more of the other columns have rational
3985
 * coefficients, but the parameter coefficients are all integral, then
3986
 * we can perform a regular (non-parametric) cut.
3987
 * Finally, if there is any parameter coefficient that is non-integral,
3988
 * then we need to involve the context tableau.  There are two cases here.
3989
 * If at least one other column has a rational coefficient, then we
3990
 * can perform a parametric cut in the main tableau by adding a new
3991
 * integer division in the context tableau.
3992
 * If all other columns have integral coefficients, then we need to
3993
 * enforce that the rational combination of parameters (c + \sum a_i y_i)/m
3994
 * is always integral.  We do this by introducing an integer division
3995
 * q = floor((c + \sum a_i y_i)/m) and stipulating that its argument should
3996
 * always be integral in the context tableau, i.e., m q = c + \sum a_i y_i.
3997
 * Since q is expressed in the tableau as
3998
 *  c + \sum a_i y_i - m q >= 0
3999
 *  -c - \sum a_i y_i + m q + m - 1 >= 0
4000
 * it is sufficient to add the inequality
4001
 *  -c - \sum a_i y_i + m q >= 0
4002
 * In the part of the context where this inequality does not hold, the
4003
 * main tableau is marked as being empty.
4004
 */
4005
static void find_solutions(struct isl_sol *sol, struct isl_tab *tab)
4006
3.95k
{
4007
3.95k
  struct isl_context *context;
4008
3.95k
  int r;
4009
3.95k
4010
3.95k
  if (
!tab || 3.95k
sol->error3.95k
)
4011
0
    goto error;
4012
3.95k
4013
3.95k
  context = sol->context;
4014
3.95k
4015
3.95k
  if (tab->empty)
4016
0
    goto done;
4017
3.95k
  
if (3.95k
context->op->is_empty(context)3.95k
)
4018
0
    goto done;
4019
3.95k
4020
5.92k
  
for (r = 0; 3.95k
r >= 0 && 5.92k
tab5.92k
&&
!tab->empty5.92k
;
r = restore_lexmin(tab)1.97k
)
{4.93k
4021
4.93k
    int flags;
4022
4.93k
    int row;
4023
4.93k
    enum isl_tab_row_sign sgn;
4024
4.93k
    int split = -1;
4025
4.93k
    int n_split = 0;
4026
4.93k
4027
39.6k
    for (row = tab->n_redundant; 
row < tab->n_row39.6k
;
++row34.7k
)
{34.8k
4028
34.8k
      if (!isl_tab_var_from_row(tab, row)->is_nonneg)
4029
1.19k
        continue;
4030
33.6k
      sgn = row_sign(tab, sol, row);
4031
33.6k
      if (!sgn)
4032
0
        goto error;
4033
33.6k
      tab->row_sign[row] = sgn;
4034
33.6k
      if (sgn == isl_tab_row_any)
4035
1.50k
        n_split++;
4036
33.6k
      if (
sgn == isl_tab_row_any && 33.6k
split == -11.50k
)
4037
1.17k
        split = row;
4038
33.6k
      if (sgn == isl_tab_row_neg)
4039
99
        break;
4040
33.6k
    }
4041
4.93k
    
if (4.93k
row < tab->n_row4.93k
)
4042
99
      continue;
4043
4.83k
    
if (4.83k
split != -14.83k
)
{1.14k
4044
1.14k
      struct isl_vec *ineq;
4045
1.14k
      if (n_split != 1)
4046
235
        split = context->op->best_split(context, tab);
4047
1.14k
      if (split < 0)
4048
0
        goto error;
4049
1.14k
      ineq = get_row_parameter_ineq(tab, split);
4050
1.14k
      if (!ineq)
4051
0
        goto error;
4052
1.14k
      is_strict(ineq);
4053
1.14k
      reset_any_to_unknown(tab);
4054
1.14k
      tab->row_sign[split] = isl_tab_row_pos;
4055
1.14k
      sol_inc_level(sol);
4056
1.14k
      find_in_pos(sol, tab, ineq->el);
4057
1.14k
      tab->row_sign[split] = isl_tab_row_neg;
4058
1.14k
      isl_seq_neg(ineq->el, ineq->el, ineq->size);
4059
1.14k
      isl_int_sub_ui(ineq->el[0], ineq->el[0], 1);
4060
1.14k
      if (!sol->error)
4061
1.14k
        context->op->add_ineq(context, ineq->el, 0, 1);
4062
1.14k
      isl_vec_free(ineq);
4063
1.14k
      if (sol->error)
4064
0
        goto error;
4065
1.14k
      continue;
4066
1.14k
    }
4067
3.68k
    
if (3.68k
tab->rational3.68k
)
4068
3
      break;
4069
3.68k
    row = first_non_integer_row(tab, &flags);
4070
3.68k
    if (row < 0)
4071
2.95k
      break;
4072
724
    
if (724
ISL_FL_ISSET724
(flags, I_PAR))
{323
4073
323
      if (
ISL_FL_ISSET323
(flags, I_VAR))
{0
4074
0
        if (isl_tab_mark_empty(tab) < 0)
4075
0
          goto error;
4076
0
        break;
4077
0
      }
4078
323
      row = add_cut(tab, row);
4079
401
    } else 
if (401
ISL_FL_ISSET401
(flags, I_VAR))
{36
4080
36
      struct isl_vec *div;
4081
36
      struct isl_vec *ineq;
4082
36
      int d;
4083
36
      div = get_row_split_div(tab, row);
4084
36
      if (!div)
4085
0
        goto error;
4086
36
      d = context->op->get_div(context, tab, div);
4087
36
      isl_vec_free(div);
4088
36
      if (d < 0)
4089
0
        goto error;
4090
36
      ineq = ineq_for_div(context->op->peek_basic_set(context), d);
4091
36
      if (!ineq)
4092
0
        goto error;
4093
36
      sol_inc_level(sol);
4094
36
      no_sol_in_strict(sol, tab, ineq);
4095
36
      isl_seq_neg(ineq->el, ineq->el, ineq->size);
4096
36
      context->op->add_ineq(context, ineq->el, 1, 1);
4097
36
      isl_vec_free(ineq);
4098
36
      if (
sol->error || 36
!context->op->is_ok(context)36
)
4099
0
        goto error;
4100
36
      tab = set_row_cst_to_div(tab, row, d);
4101
36
      if (context->op->is_empty(context))
4102
0
        break;
4103
36
    } else
4104
365
      row = add_parametric_cut(tab, row, context);
4105
724
    
if (724
row < 0724
)
4106
0
      goto error;
4107
724
  }
4108
3.95k
  
if (3.95k
r < 03.95k
)
4109
0
    goto error;
4110
3.95k
done:
4111
3.95k
  sol_add(sol, tab);
4112
3.95k
  isl_tab_free(tab);
4113
3.95k
  return;
4114
0
error:
4115
0
  isl_tab_free(tab);
4116
0
  sol->error = 1;
4117
0
}
4118
4119
/* Does "sol" contain a pair of partial solutions that could potentially
4120
 * be merged?
4121
 *
4122
 * We currently only check that "sol" is not in an error state
4123
 * and that there are at least two partial solutions of which the final two
4124
 * are defined at the same level.
4125
 */
4126
static int sol_has_mergeable_solutions(struct isl_sol *sol)
4127
2.80k
{
4128
2.80k
  if (sol->error)
4129
0
    return 0;
4130
2.80k
  
if (2.80k
!sol->partial2.80k
)
4131
58
    return 0;
4132
2.74k
  
if (2.74k
!sol->partial->next2.74k
)
4133
2.09k
    return 0;
4134
649
  return sol->partial->level == sol->partial->next->level;
4135
2.74k
}
4136
4137
/* Compute the lexicographic minimum of the set represented by the main
4138
 * tableau "tab" within the context "sol->context_tab".
4139
 *
4140
 * As a preprocessing step, we first transfer all the purely parametric
4141
 * equalities from the main tableau to the context tableau, i.e.,
4142
 * parameters that have been pivoted to a row.
4143
 * These equalities are ignored by the main algorithm, because the
4144
 * corresponding rows may not be marked as being non-negative.
4145
 * In parts of the context where the added equality does not hold,
4146
 * the main tableau is marked as being empty.
4147
 *
4148
 * Before we embark on the actual computation, we save a copy
4149
 * of the context.  When we return, we check if there are any
4150
 * partial solutions that can potentially be merged.  If so,
4151
 * we perform a rollback to the initial state of the context.
4152
 * The merging of partial solutions happens inside calls to
4153
 * sol_dec_level that are pushed onto the undo stack of the context.
4154
 * If there are no partial solutions that can potentially be merged
4155
 * then the rollback is skipped as it would just be wasted effort.
4156
 */
4157
static void find_solutions_main(struct isl_sol *sol, struct isl_tab *tab)
4158
2.80k
{
4159
2.80k
  int row;
4160
2.80k
  void *saved;
4161
2.80k
4162
2.80k
  if (!tab)
4163
0
    goto error;
4164
2.80k
4165
2.80k
  sol->level = 0;
4166
2.80k
4167
29.4k
  for (row = tab->n_redundant; 
row < tab->n_row29.4k
;
++row26.6k
)
{26.6k
4168
26.6k
    int p;
4169
26.6k
    struct isl_vec *eq;
4170
26.6k
4171
26.6k
    if (tab->row_var[row] < 0)
4172
4.77k
      continue;
4173
21.8k
    
if (21.8k
tab->row_var[row] >= tab->n_param &&21.8k
4174
18.3k
        tab->row_var[row] < tab->n_var - tab->n_div)
4175
18.3k
      continue;
4176
3.54k
    
if (3.54k
tab->row_var[row] < tab->n_param3.54k
)
4177
3.54k
      p = tab->row_var[row];
4178
3.54k
    else
4179
0
      p = tab->row_var[row]
4180
0
        + tab->n_param - (tab->n_var - tab->n_div);
4181
3.54k
4182
3.54k
    eq = isl_vec_alloc(tab->mat->ctx, 1+tab->n_param+tab->n_div);
4183
3.54k
    if (!eq)
4184
0
      goto error;
4185
3.54k
    get_row_parameter_line(tab, row, eq->el);
4186
3.54k
    isl_int_neg(eq->el[1 + p], tab->mat->row[row][0]);
4187
3.54k
    eq = isl_vec_normalize(eq);
4188
3.54k
4189
3.54k
    sol_inc_level(sol);
4190
3.54k
    no_sol_in_strict(sol, tab, eq);
4191
3.54k
4192
3.54k
    isl_seq_neg(eq->el, eq->el, eq->size);
4193
3.54k
    sol_inc_level(sol);
4194
3.54k
    no_sol_in_strict(sol, tab, eq);
4195
3.54k
    isl_seq_neg(eq->el, eq->el, eq->size);
4196
3.54k
4197
3.54k
    sol->context->op->add_eq(sol->context, eq->el, 1, 1);
4198
3.54k
4199
3.54k
    isl_vec_free(eq);
4200
3.54k
4201
3.54k
    if (isl_tab_mark_redundant(tab, row) < 0)
4202
0
      goto error;
4203
3.54k
4204
3.54k
    
if (3.54k
sol->context->op->is_empty(sol->context)3.54k
)
4205
0
      break;
4206
3.54k
4207
3.54k
    row = tab->n_redundant - 1;
4208
3.54k
  }
4209
2.80k
4210
2.80k
  saved = sol->context->op->save(sol->context);
4211
2.80k
4212
2.80k
  find_solutions(sol, tab);
4213
2.80k
4214
2.80k
  if (sol_has_mergeable_solutions(sol))
4215
620
    sol->context->op->restore(sol->context, saved);
4216
2.80k
  else
4217
2.18k
    sol->context->op->discard(saved);
4218
2.80k
4219
2.80k
  sol->level = 0;
4220
2.80k
  sol_pop(sol);
4221
2.80k
4222
2.80k
  return;
4223
0
error:
4224
0
  isl_tab_free(tab);
4225
0
  sol->error = 1;
4226
0
}
4227
4228
/* Check if integer division "div" of "dom" also occurs in "bmap".
4229
 * If so, return its position within the divs.
4230
 * If not, return -1.
4231
 */
4232
static int find_context_div(struct isl_basic_map *bmap,
4233
  struct isl_basic_set *dom, unsigned div)
4234
204
{
4235
204
  int i;
4236
204
  unsigned b_dim = isl_space_dim(bmap->dim, isl_dim_all);
4237
204
  unsigned d_dim = isl_space_dim(dom->dim, isl_dim_all);
4238
204
4239
204
  if (isl_int_is_zero(dom->div[div][0]))
4240
10
    return -1;
4241
194
  
if (194
isl_seq_first_non_zero(dom->div[div] + 2 + d_dim, dom->n_div) != -1194
)
4242
2
    return -1;
4243
194
4244
242
  
for (i = 0; 192
i < bmap->n_div242
;
++i50
)
{74
4245
74
    if (isl_int_is_zero(bmap->div[i][0]))
4246
10
      continue;
4247
64
    
if (64
isl_seq_first_non_zero(bmap->div[i] + 2 + d_dim,64
4248
64
             (b_dim - d_dim) + bmap->n_div) != -1)
4249
10
      continue;
4250
54
    
if (54
isl_seq_eq(bmap->div[i], dom->div[div], 2 + d_dim)54
)
4251
24
      return i;
4252
54
  }
4253
168
  return -1;
4254
192
}
4255
4256
/* The correspondence between the variables in the main tableau,
4257
 * the context tableau, and the input map and domain is as follows.
4258
 * The first n_param and the last n_div variables of the main tableau
4259
 * form the variables of the context tableau.
4260
 * In the basic map, these n_param variables correspond to the
4261
 * parameters and the input dimensions.  In the domain, they correspond
4262
 * to the parameters and the set dimensions.
4263
 * The n_div variables correspond to the integer divisions in the domain.
4264
 * To ensure that everything lines up, we may need to copy some of the
4265
 * integer divisions of the domain to the map.  These have to be placed
4266
 * in the same order as those in the context and they have to be placed
4267
 * after any other integer divisions that the map may have.
4268
 * This function performs the required reordering.
4269
 */
4270
static struct isl_basic_map *align_context_divs(struct isl_basic_map *bmap,
4271
  struct isl_basic_set *dom)
4272
95
{
4273
95
  int i;
4274
95
  int common = 0;
4275
95
  int other;
4276
95
4277
197
  for (i = 0; 
i < dom->n_div197
;
++i102
)
4278
102
    
if (102
find_context_div(bmap, dom, i) != -1102
)
4279
12
      common++;
4280
95
  other = bmap->n_div - common;
4281
95
  if (
dom->n_div - common > 095
)
{85
4282
85
    bmap = isl_basic_map_extend_space(bmap, isl_space_copy(bmap->dim),
4283
85
        dom->n_div - common, 0, 0);
4284
85
    if (!bmap)
4285
0
      return NULL;
4286
85
  }
4287
197
  
for (i = 0; 95
i < dom->n_div197
;
++i102
)
{102
4288
102
    int pos = find_context_div(bmap, dom, i);
4289
102
    if (
pos < 0102
)
{90
4290
90
      pos = isl_basic_map_alloc_div(bmap);
4291
90
      if (pos < 0)
4292
0
        goto error;
4293
90
      
isl_int_set_si90
(bmap->div[pos][0], 0);90
4294
90
    }
4295
102
    
if (102
pos != other + i102
)
4296
7
      isl_basic_map_swap_div(bmap, pos, other + i);
4297
102
  }
4298
95
  return bmap;
4299
0
error:
4300
0
  isl_basic_map_free(bmap);
4301
0
  return NULL;
4302
95
}
4303
4304
/* Base case of isl_tab_basic_map_partial_lexopt, after removing
4305
 * some obvious symmetries.
4306
 *
4307
 * We make sure the divs in the domain are properly ordered,
4308
 * because they will be added one by one in the given order
4309
 * during the construction of the solution map.
4310
 * Furthermore, make sure that the known integer divisions
4311
 * appear before any unknown integer division because the solution
4312
 * may depend on the known integer divisions, while anything that
4313
 * depends on any variable starting from the first unknown integer
4314
 * division is ignored in sol_pma_add.
4315
 */
4316
static struct isl_sol *basic_map_partial_lexopt_base_sol(
4317
  __isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom,
4318
  __isl_give isl_set **empty, int max,
4319
  struct isl_sol *(*init)(__isl_keep isl_basic_map *bmap,
4320
        __isl_take isl_basic_set *dom, int track_empty, int max))
4321
3.04k
{
4322
3.04k
  struct isl_tab *tab;
4323
3.04k
  struct isl_sol *sol = NULL;
4324
3.04k
  struct isl_context *context;
4325
3.04k
4326
3.04k
  if (
dom->n_div3.04k
)
{95
4327
95
    dom = isl_basic_set_sort_divs(dom);
4328
95
    bmap = align_context_divs(bmap, dom);
4329
95
  }
4330
3.04k
  sol = init(bmap, dom, !!empty, max);
4331
3.04k
  if (!sol)
4332
0
    goto error;
4333
3.04k
4334
3.04k
  context = sol->context;
4335
3.04k
  if (isl_basic_set_plain_is_empty(context->op->peek_basic_set(context)))
4336
0
    /* nothing */;
4337
3.04k
  else 
if (3.04k
isl_basic_map_plain_is_empty(bmap)3.04k
)
{420
4338
420
    if (sol->add_empty)
4339
420
      sol->add_empty(sol,
4340
420
        isl_basic_set_copy(context->op->peek_basic_set(context)));
4341
2.62k
  } else {
4342
2.62k
    tab = tab_for_lexmin(bmap,
4343
2.62k
            context->op->peek_basic_set(context), 1, max);
4344
2.62k
    tab = context->op->detect_nonnegative_parameters(context, tab);
4345
2.62k
    find_solutions_main(sol, tab);
4346
2.62k
  }
4347
3.04k
  if (sol->error)
4348
0
    goto error;
4349
3.04k
4350
3.04k
  isl_basic_map_free(bmap);
4351
3.04k
  return sol;
4352
0
error:
4353
0
  sol_free(sol);
4354
0
  isl_basic_map_free(bmap);
4355
0
  return NULL;
4356
3.04k
}
4357
4358
/* Base case of isl_tab_basic_map_partial_lexopt, after removing
4359
 * some obvious symmetries.
4360
 *
4361
 * We call basic_map_partial_lexopt_base_sol and extract the results.
4362
 */
4363
static __isl_give isl_map *basic_map_partial_lexopt_base(
4364
  __isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom,
4365
  __isl_give isl_set **empty, int max)
4366
1.27k
{
4367
1.27k
  isl_map *result = NULL;
4368
1.27k
  struct isl_sol *sol;
4369
1.27k
  struct isl_sol_map *sol_map;
4370
1.27k
4371
1.27k
  sol = basic_map_partial_lexopt_base_sol(bmap, dom, empty, max,
4372
1.27k
            &sol_map_init);
4373
1.27k
  if (!sol)
4374
0
    return NULL;
4375
1.27k
  sol_map = (struct isl_sol_map *) sol;
4376
1.27k
4377
1.27k
  result = isl_map_copy(sol_map->map);
4378
1.27k
  if (empty)
4379
941
    *empty = isl_set_copy(sol_map->empty);
4380
1.27k
  sol_free(&sol_map->sol);
4381
1.27k
  return result;
4382
1.27k
}
4383
4384
/* Return a count of the number of occurrences of the "n" first
4385
 * variables in the inequality constraints of "bmap".
4386
 */
4387
static __isl_give int *count_occurrences(__isl_keep isl_basic_map *bmap,
4388
  int n)
4389
3.06k
{
4390
3.06k
  int i, j;
4391
3.06k
  isl_ctx *ctx;
4392
3.06k
  int *occurrences;
4393
3.06k
4394
3.06k
  if (!bmap)
4395
0
    return NULL;
4396
3.06k
  ctx = isl_basic_map_get_ctx(bmap);
4397
3.06k
  occurrences = isl_calloc_array(ctx, int, n);
4398
3.06k
  if (!occurrences)
4399
0
    return NULL;
4400
3.06k
4401
9.36k
  
for (i = 0; 3.06k
i < bmap->n_ineq9.36k
;
++i6.29k
)
{6.29k
4402
33.1k
    for (j = 0; 
j < n33.1k
;
++j26.8k
)
{26.8k
4403
26.8k
      if (
!26.8k
isl_int_is_zero26.8k
(bmap->ineq[i][1 + j]))
4404
4.91k
        occurrences[j]++;
4405
26.8k
    }
4406
6.29k
  }
4407
3.06k
4408
3.06k
  return occurrences;
4409
3.06k
}
4410
4411
/* Do all of the "n" variables with non-zero coefficients in "c"
4412
 * occur in exactly a single constraint.
4413
 * "occurrences" is an array of length "n" containing the number
4414
 * of occurrences of each of the variables in the inequality constraints.
4415
 */
4416
static int single_occurrence(int n, isl_int *c, int *occurrences)
4417
3.13k
{
4418
3.13k
  int i;
4419
3.13k
4420
10.5k
  for (i = 0; 
i < n10.5k
;
++i7.46k
)
{8.17k
4421
8.17k
    if (isl_int_is_zero(c[i]))
4422
7.15k
      continue;
4423
1.02k
    
if (1.02k
occurrences[i] != 11.02k
)
4424
710
      return 0;
4425
1.02k
  }
4426
3.13k
4427
2.42k
  return 1;
4428
3.13k
}
4429
4430
/* Do all of the "n" initial variables that occur in inequality constraint
4431
 * "ineq" of "bmap" only occur in that constraint?
4432
 */
4433
static int all_single_occurrence(__isl_keep isl_basic_map *bmap, int ineq,
4434
  int n)
4435
0
{
4436
0
  int i, j;
4437
0
4438
0
  for (i = 0; 
i < n0
;
++i0
)
{0
4439
0
    if (isl_int_is_zero(bmap->ineq[ineq][1 + i]))
4440
0
      continue;
4441
0
    
for (j = 0; 0
j < bmap->n_ineq0
;
++j0
)
{0
4442
0
      if (j == ineq)
4443
0
        continue;
4444
0
      
if (0
!0
isl_int_is_zero0
(bmap->ineq[j][1 + i]))
4445
0
        return 0;
4446
0
    }
4447
0
  }
4448
0
4449
0
  return 1;
4450
0
}
4451
4452
/* Structure used during detection of parallel constraints.
4453
 * n_in: number of "input" variables: isl_dim_param + isl_dim_in
4454
 * n_out: number of "output" variables: isl_dim_out + isl_dim_div
4455
 * val: the coefficients of the output variables
4456
 */
4457
struct isl_constraint_equal_info {
4458
  unsigned n_in;
4459
  unsigned n_out;
4460
  isl_int *val;
4461
};
4462
4463
/* Check whether the coefficients of the output variables
4464
 * of the constraint in "entry" are equal to info->val.
4465
 */
4466
static int constraint_equal(const void *entry, const void *val)
4467
27
{
4468
27
  isl_int **row = (isl_int **)entry;
4469
27
  const struct isl_constraint_equal_info *info = val;
4470
27
4471
27
  return isl_seq_eq((*row) + 1 + info->n_in, info->val, info->n_out);
4472
27
}
4473
4474
/* Check whether "bmap" has a pair of constraints that have
4475
 * the same coefficients for the output variables.
4476
 * Note that the coefficients of the existentially quantified
4477
 * variables need to be zero since the existentially quantified
4478
 * of the result are usually not the same as those of the input.
4479
 * Furthermore, check that each of the input variables that occur
4480
 * in those constraints does not occur in any other constraint.
4481
 * If so, return true and return the row indices of the two constraints
4482
 * in *first and *second.
4483
 */
4484
static isl_bool parallel_constraints(__isl_keep isl_basic_map *bmap,
4485
  int *first, int *second)
4486
3.06k
{
4487
3.06k
  int i;
4488
3.06k
  isl_ctx *ctx;
4489
3.06k
  int *occurrences = NULL;
4490
3.06k
  struct isl_hash_table *table = NULL;
4491
3.06k
  struct isl_hash_table_entry *entry;
4492
3.06k
  struct isl_constraint_equal_info info;
4493
3.06k
  unsigned n_out;
4494
3.06k
  unsigned n_div;
4495
3.06k
4496
3.06k
  ctx = isl_basic_map_get_ctx(bmap);
4497
3.06k
  table = isl_hash_table_alloc(ctx, bmap->n_ineq);
4498
3.06k
  if (!table)
4499
0
    goto error;
4500
3.06k
4501
3.06k
  info.n_in = isl_basic_map_dim(bmap, isl_dim_param) +
4502
3.06k
        isl_basic_map_dim(bmap, isl_dim_in);
4503
3.06k
  occurrences = count_occurrences(bmap, info.n_in);
4504
3.06k
  if (
info.n_in && 3.06k
!occurrences2.61k
)
4505
0
    goto error;
4506
3.06k
  n_out = isl_basic_map_dim(bmap, isl_dim_out);
4507
3.06k
  n_div = isl_basic_map_dim(bmap, isl_dim_div);
4508
3.06k
  info.n_out = n_out + n_div;
4509
9.30k
  for (i = 0; 
i < bmap->n_ineq9.30k
;
++i6.23k
)
{6.26k
4510
6.26k
    uint32_t hash;
4511
6.26k
4512
6.26k
    info.val = bmap->ineq[i] + 1 + info.n_in;
4513
6.26k
    if (isl_seq_first_non_zero(info.val, n_out) < 0)
4514
3.10k
      continue;
4515
3.16k
    
if (3.16k
isl_seq_first_non_zero(info.val + n_out, n_div) >= 03.16k
)
4516
32
      continue;
4517
3.13k
    
if (3.13k
!single_occurrence(info.n_in, bmap->ineq[i] + 1,3.13k
4518
3.13k
          occurrences))
4519
710
      continue;
4520
2.42k
    hash = isl_seq_get_hash(info.val, info.n_out);
4521
2.42k
    entry = isl_hash_table_find(ctx, table, hash,
4522
2.42k
              constraint_equal, &info, 1);
4523
2.42k
    if (!entry)
4524
0
      goto error;
4525
2.42k
    
if (2.42k
entry->data2.42k
)
4526
27
      break;
4527
2.39k
    entry->data = &bmap->ineq[i];
4528
2.39k
  }
4529
3.06k
4530
3.06k
  
if (3.06k
i < bmap->n_ineq3.06k
)
{27
4531
27
    *first = ((isl_int **)entry->data) - bmap->ineq; 
4532
27
    *second = i;
4533
27
  }
4534
3.06k
4535
3.06k
  isl_hash_table_free(ctx, table);
4536
3.06k
  free(occurrences);
4537
3.06k
4538
3.06k
  return i < bmap->n_ineq;
4539
0
error:
4540
0
  isl_hash_table_free(ctx, table);
4541
0
  free(occurrences);
4542
0
  return isl_bool_error;
4543
3.06k
}
4544
4545
/* Given a set of upper bounds in "var", add constraints to "bset"
4546
 * that make the i-th bound smallest.
4547
 *
4548
 * In particular, if there are n bounds b_i, then add the constraints
4549
 *
4550
 *  b_i <= b_j  for j > i
4551
 *  b_i <  b_j  for j < i
4552
 */
4553
static __isl_give isl_basic_set *select_minimum(__isl_take isl_basic_set *bset,
4554
  __isl_keep isl_mat *var, int i)
4555
90
{
4556
90
  isl_ctx *ctx;
4557
90
  int j, k;
4558
90
4559
90
  ctx = isl_mat_get_ctx(var);
4560
90
4561
270
  for (j = 0; 
j < var->n_row270
;
++j180
)
{180
4562
180
    if (j == i)
4563
90
      continue;
4564
90
    k = isl_basic_set_alloc_inequality(bset);
4565
90
    if (k < 0)
4566
0
      goto error;
4567
90
    isl_seq_combine(bset->ineq[k], ctx->one, var->row[j],
4568
90
        ctx->negone, var->row[i], var->n_col);
4569
90
    isl_int_set_si(bset->ineq[k][var->n_col], 0);
4570
90
    if (j < i)
4571
45
      isl_int_sub_ui(bset->ineq[k][0], bset->ineq[k][0], 1);
4572
90
  }
4573
90
4574
90
  bset = isl_basic_set_finalize(bset);
4575
90
4576
90
  return bset;
4577
0
error:
4578
0
  isl_basic_set_free(bset);
4579
0
  return NULL;
4580
90
}
4581
4582
/* Given a set of upper bounds on the last "input" variable m,
4583
 * construct a set that assigns the minimal upper bound to m, i.e.,
4584
 * construct a set that divides the space into cells where one
4585
 * of the upper bounds is smaller than all the others and assign
4586
 * this upper bound to m.
4587
 *
4588
 * In particular, if there are n bounds b_i, then the result
4589
 * consists of n basic sets, each one of the form
4590
 *
4591
 *  m = b_i
4592
 *  b_i <= b_j  for j > i
4593
 *  b_i <  b_j  for j < i
4594
 */
4595
static __isl_give isl_set *set_minimum(__isl_take isl_space *dim,
4596
  __isl_take isl_mat *var)
4597
27
{
4598
27
  int i, k;
4599
27
  isl_basic_set *bset = NULL;
4600
27
  isl_set *set = NULL;
4601
27
4602
27
  if (
!dim || 27
!var27
)
4603
0
    goto error;
4604
27
4605
27
  set = isl_set_alloc_space(isl_space_copy(dim),
4606
27
        var->n_row, ISL_SET_DISJOINT);
4607
27
4608
81
  for (i = 0; 
i < var->n_row81
;
++i54
)
{54
4609
54
    bset = isl_basic_set_alloc_space(isl_space_copy(dim), 0,
4610
54
                 1, var->n_row - 1);
4611
54
    k = isl_basic_set_alloc_equality(bset);
4612
54
    if (k < 0)
4613
0
      goto error;
4614
54
    isl_seq_cpy(bset->eq[k], var->row[i], var->n_col);
4615
54
    isl_int_set_si(bset->eq[k][var->n_col], -1);
4616
54
    bset = select_minimum(bset, var, i);
4617
54
    set = isl_set_add_basic_set(set, bset);
4618
54
  }
4619
27
4620
27
  isl_space_free(dim);
4621
27
  isl_mat_free(var);
4622
27
  return set;
4623
0
error:
4624
0
  isl_basic_set_free(bset);
4625
0
  isl_set_free(set);
4626
0
  isl_space_free(dim);
4627
0
  isl_mat_free(var);
4628
0
  return NULL;
4629
27
}
4630
4631
/* Given that the last input variable of "bmap" represents the minimum
4632
 * of the bounds in "cst", check whether we need to split the domain
4633
 * based on which bound attains the minimum.
4634
 *
4635
 * A split is needed when the minimum appears in an integer division
4636
 * or in an equality.  Otherwise, it is only needed if it appears in
4637
 * an upper bound that is different from the upper bounds on which it
4638
 * is defined.
4639
 */
4640
static isl_bool need_split_basic_map(__isl_keep isl_basic_map *bmap,
4641
  __isl_keep isl_mat *cst)
4642
13
{
4643
13
  int i, j;
4644
13
  unsigned total;
4645
13
  unsigned pos;
4646
13
4647
13
  pos = cst->n_col - 1;
4648
13
  total = isl_basic_map_dim(bmap, isl_dim_all);
4649
13
4650
17
  for (i = 0; 
i < bmap->n_div17
;
++i4
)
4651
8
    
if (8
!8
isl_int_is_zero8
(bmap->div[i][2 + pos]))
4652
4
      return isl_bool_true;
4653
13
4654
21
  
for (i = 0; 9
i < bmap->n_eq21
;
++i12
)
4655
13
    
if (13
!13
isl_int_is_zero13
(bmap->eq[i][1 + pos]))
4656
1
      return isl_bool_true;
4657
9
4658
30
  
for (i = 0; 8
i < bmap->n_ineq30
;
++i22
)
{28
4659
28
    if (isl_int_is_nonneg(bmap->ineq[i][1 + pos]))
4660
14
      continue;
4661
14
    
if (14
!14
isl_int_is_negone14
(bmap->ineq[i][1 + pos]))
4662
0
      return isl_bool_true;
4663
14
    
if (14
isl_seq_first_non_zero(bmap->ineq[i] + 1 + pos + 1,14
4664
14
             total - pos - 1) >= 0)
4665
4
      return isl_bool_true;
4666
14
4667
18
    
for (j = 0; 10
j < cst->n_row18
;
++j8
)
4668
16
      
if (16
isl_seq_eq(bmap->ineq[i], cst->row[j], cst->n_col)16
)
4669
8
        break;
4670
10
    if (j >= cst->n_row)
4671
2
      return isl_bool_true;
4672
10
  }
4673
8
4674
2
  return isl_bool_false;
4675
8
}
4676
4677
/* Given that the last set variable of "bset" represents the minimum
4678
 * of the bounds in "cst", check whether we need to split the domain
4679
 * based on which bound attains the minimum.
4680
 *
4681
 * We simply call need_split_basic_map here.  This is safe because
4682
 * the position of the minimum is computed from "cst" and not
4683
 * from "bmap".
4684
 */
4685
static isl_bool need_split_basic_set(__isl_keep isl_basic_set *bset,
4686
  __isl_keep isl_mat *cst)
4687
2
{
4688
2
  return need_split_basic_map(bset_to_bmap(bset), cst);
4689
2
}
4690
4691
/* Given that the last set variable of "set" represents the minimum
4692
 * of the bounds in "cst", check whether we need to split the domain
4693
 * based on which bound attains the minimum.
4694
 */
4695
static isl_bool need_split_set(__isl_keep isl_set *set, __isl_keep isl_mat *cst)
4696
0
{
4697
0
  int i;
4698
0
4699
0
  for (i = 0; 
i < set->n0
;
++i0
)
{0
4700
0
    isl_bool split;
4701
0
4702
0
    split = need_split_basic_set(set->p[i], cst);
4703
0
    if (
split < 0 || 0
split0
)
4704
0
      return split;
4705
0
  }
4706
0
4707
0
  return isl_bool_false;
4708
0
}
4709
4710
/* Given a set of which the last set variable is the minimum
4711
 * of the bounds in "cst", split each basic set in the set
4712
 * in pieces where one of the bounds is (strictly) smaller than the others.
4713
 * This subdivision is given in "min_expr".
4714
 * The variable is subsequently projected out.
4715
 *
4716
 * We only do the split when it is needed.
4717
 * For example if the last input variable m = min(a,b) and the only
4718
 * constraints in the given basic set are lower bounds on m,
4719
 * i.e., l <= m = min(a,b), then we can simply project out m
4720
 * to obtain l <= a and l <= b, without having to split on whether
4721
 * m is equal to a or b.
4722
 */
4723
static __isl_give isl_set *split(__isl_take isl_set *empty,
4724
  __isl_take isl_set *min_expr, __isl_take isl_mat *cst)
4725
2
{
4726
2
  int n_in;
4727
2
  int i;
4728
2
  isl_space *dim;
4729
2
  isl_set *res;
4730
2
4731
2
  if (
!empty || 2
!min_expr2
||
!cst2
)
4732
0
    goto error;
4733
2
4734
2
  n_in = isl_set_dim(empty, isl_dim_set);
4735
2
  dim = isl_set_get_space(empty);
4736
2
  dim = isl_space_drop_dims(dim, isl_dim_set, n_in - 1, 1);
4737
2
  res = isl_set_empty(dim);
4738
2
4739
4
  for (i = 0; 
i < empty->n4
;
++i2
)
{2
4740
2
    isl_bool split;
4741
2
    isl_set *set;
4742
2
4743
2
    set = isl_set_from_basic_set(isl_basic_set_copy(empty->p[i]));
4744
2
    split = need_split_basic_set(empty->p[i], cst);
4745
2
    if (split < 0)
4746
0
      set = isl_set_free(set);
4747
2
    else 
if (2
split2
)
4748
2
      set = isl_set_intersect(set, isl_set_copy(min_expr));
4749
2
    set = isl_set_remove_dims(set, isl_dim_set, n_in - 1, 1);
4750
2
4751
2
    res = isl_set_union_disjoint(res, set);
4752
2
  }
4753
2
4754
2
  isl_set_free(empty);
4755
2
  isl_set_free(min_expr);
4756
2
  isl_mat_free(cst);
4757
2
  return res;
4758
0
error:
4759
0
  isl_set_free(empty);
4760
0
  isl_set_free(min_expr);
4761
0
  isl_mat_free(cst);
4762
0
  return NULL;
4763
2
}
4764
4765
/* Given a map of which the last input variable is the minimum
4766
 * of the bounds in "cst", split each basic set in the set
4767
 * in pieces where one of the bounds is (strictly) smaller than the others.
4768
 * This subdivision is given in "min_expr".
4769
 * The variable is subsequently projected out.
4770
 *
4771
 * The implementation is essentially the same as that of "split".
4772
 */
4773
static __isl_give isl_map *split_domain(__isl_take isl_map *opt,
4774
  __isl_take isl_set *min_expr, __isl_take isl_mat *cst)
4775
9
{
4776
9
  int n_in;
4777
9
  int i;
4778
9
  isl_space *dim;
4779
9
  isl_map *res;
4780
9
4781
9
  if (
!opt || 9
!min_expr9
||
!cst9
)
4782
0
    goto error;
4783
9
4784
9
  n_in = isl_map_dim(opt, isl_dim_in);
4785
9
  dim = isl_map_get_space(opt);
4786
9
  dim = isl_space_drop_dims(dim, isl_dim_in, n_in - 1, 1);
4787
9
  res = isl_map_empty(dim);
4788
9
4789
20
  for (i = 0; 
i < opt->n20
;
++i11
)
{11
4790
11
    isl_map *map;
4791
11
    isl_bool split;
4792
11
4793
11
    map = isl_map_from_basic_map(isl_basic_map_copy(opt->p[i]));
4794
11
    split = need_split_basic_map(opt->p