Coverage Report

Created: 2017-04-29 12:21

/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
4.76k
{
202
4.76k
  struct isl_partial_sol *partial, *next;
203
4.76k
  if (!sol)
204
0
    return;
205
4.76k
  
for (partial = sol->partial; 4.76k
partial4.76k
;
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
4.76k
  isl_space_free(sol->space);
212
4.76k
  if (sol->context)
213
4.76k
    sol->context->op->free(sol->context);
214
4.76k
  sol->free(sol);
215
4.76k
  free(sol);
216
4.76k
}
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
5.79k
{
226
5.79k
  struct isl_partial_sol *partial;
227
5.79k
228
5.79k
  if (
sol->error || 5.79k
!dom5.79k
)
229
0
    goto error;
230
5.79k
231
5.79k
  
partial = 5.79k
isl_alloc_type5.79k
(dom->ctx, struct isl_partial_sol);
232
5.79k
  if (!partial)
233
0
    goto error;
234
5.79k
235
5.79k
  partial->level = sol->level;
236
5.79k
  partial->dom = dom;
237
5.79k
  partial->ma = ma;
238
5.79k
  partial->next = sol->partial;
239
5.79k
240
5.79k
  sol->partial = partial;
241
5.79k
242
5.79k
  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
3.96k
{
254
3.96k
  int i;
255
3.96k
  unsigned rows, cols, n;
256
3.96k
257
3.96k
  if (!M)
258
0
    return isl_stat_error;
259
3.96k
  rows = isl_mat_rows(M);
260
3.96k
  cols = isl_mat_cols(M);
261
3.96k
  n = cols - first;
262
17.5k
  for (i = 0; 
i < rows17.5k
;
++i13.5k
)
263
13.5k
    
if (13.5k
isl_seq_first_non_zero(M->row[i] + first, n) != -113.5k
)
264
0
      isl_die(isl_mat_get_ctx(M), isl_error_internal,
265
3.96k
        "final columns should be zero",
266
3.96k
        return isl_stat_error);
267
3.96k
  return isl_stat_ok;
268
3.96k
}
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
3.96k
{
279
3.96k
  int i, dim;
280
3.96k
  isl_aff *aff;
281
3.96k
282
3.96k
  if (
!ma || 3.96k
!ls3.96k
||
!M3.96k
)
283
0
    goto error;
284
3.96k
285
3.96k
  dim = isl_local_space_dim(ls, isl_dim_all);
286
3.96k
  if (check_final_columns_are_zero(M, 1 + dim) < 0)
287
0
    goto error;
288
13.5k
  
for (i = 1; 3.96k
i < M->n_row13.5k
;
++i9.57k
)
{9.57k
289
9.57k
    aff = isl_aff_alloc(isl_local_space_copy(ls));
290
9.57k
    if (
aff9.57k
)
{9.57k
291
9.57k
      isl_int_set(aff->v->el[0], M->row[0][0]);
292
9.57k
      isl_seq_cpy(aff->v->el + 1, M->row[i], 1 + dim);
293
9.57k
    }
294
9.57k
    aff = isl_aff_normalize(aff);
295
9.57k
    ma = isl_multi_aff_set_aff(ma, i - 1, aff);
296
9.57k
  }
297
3.96k
  isl_local_space_free(ls);
298
3.96k
  isl_mat_free(M);
299
3.96k
300
3.96k
  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
3.96k
}
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
3.96k
{
331
3.96k
  isl_local_space *ls;
332
3.96k
  isl_multi_aff *ma;
333
3.96k
  int n_div, n_known;
334
3.96k
335
3.96k
  n_div = isl_basic_set_dim(dom, isl_dim_div);
336
3.96k
  n_known = n_div - sol->context->n_unknown;
337
3.96k
338
3.96k
  ma = isl_multi_aff_alloc(isl_space_copy(sol->space));
339
3.96k
  ls = isl_basic_set_get_local_space(dom);
340
3.96k
  ls = isl_local_space_drop_dims(ls, isl_dim_div,
341
3.96k
          n_known, n_div - n_known);
342
3.96k
  ma = set_from_affine_matrix(ma, ls, M);
343
3.96k
344
3.96k
  if (!ma)
345
0
    dom = isl_basic_set_free(dom);
346
3.96k
  sol_push_sol(sol, dom, ma);
347
3.96k
}
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
5.71k
{
354
5.71k
  struct isl_partial_sol *partial;
355
5.71k
356
5.71k
  partial = sol->partial;
357
5.71k
  sol->partial = partial->next;
358
5.71k
359
5.71k
  if (partial->ma)
360
3.88k
    sol->add(sol, partial->dom, partial->ma);
361
5.71k
  else
362
1.82k
    sol->add_empty(sol, partial->dom);
363
5.71k
  free(partial);
364
5.71k
}
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
5.88k
{
370
5.88k
  struct isl_basic_set *bset;
371
5.88k
372
5.88k
  if (sol->error)
373
0
    return NULL;
374
5.88k
375
5.88k
  bset = isl_basic_set_dup(sol->context->op->peek_basic_set(sol->context));
376
5.88k
  bset = isl_basic_set_update_from_tab(bset,
377
5.88k
      sol->context->op->peek_tab(sol->context));
378
5.88k
379
5.88k
  return bset;
380
5.88k
}
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
1.48k
{
387
1.48k
  if (!s1->ma != !s2->ma)
388
1.38k
    return isl_bool_false;
389
98
  
if (98
!s1->ma98
)
390
0
    return isl_bool_true;
391
98
392
98
  return isl_multi_aff_plain_is_equal(s1->ma, s2->ma);
393
98
}
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
84
{
427
84
  struct isl_partial_sol *partial;
428
84
  isl_basic_set *bset;
429
84
430
84
  partial = sol->partial;
431
84
432
84
  bset = sol_domain(sol);
433
84
  isl_basic_set_free(partial->next->dom);
434
84
  partial->next->dom = bset;
435
84
  partial->next->level = sol->level;
436
84
437
84
  if (!bset)
438
0
    return isl_stat_error;
439
84
440
84
  sol->partial = partial->next;
441
84
  isl_basic_set_free(partial->dom);
442
84
  isl_multi_aff_free(partial->ma);
443
84
  free(partial);
444
84
445
84
  return isl_stat_ok;
446
84
}
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
139
{
458
139
  isl_set *set;
459
139
  isl_pw_multi_aff *pma1, *pma2;
460
139
  isl_bool equal;
461
139
462
139
  set = isl_basic_set_compute_divs(isl_basic_set_copy(dom));
463
139
  pma1 = isl_pw_multi_aff_alloc(isl_set_copy(set),
464
139
          isl_multi_aff_copy(ma1));
465
139
  pma2 = isl_pw_multi_aff_alloc(set, isl_multi_aff_copy(ma2));
466
139
  equal = isl_pw_multi_aff_is_equal(pma1, pma2);
467
139
  isl_pw_multi_aff_free(pma1);
468
139
  isl_pw_multi_aff_free(pma2);
469
139
470
139
  return equal;
471
139
}
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
1.48k
{
490
1.48k
  struct isl_partial_sol *partial;
491
1.48k
  isl_bool same;
492
1.48k
493
1.48k
  partial = sol->partial;
494
1.48k
495
1.48k
  same = same_solution(partial, partial->next);
496
1.48k
  if (same < 0)
497
0
    return isl_stat_error;
498
1.48k
  
if (1.48k
same1.48k
)
499
9
    return combine_initial_into_second(sol);
500
1.47k
  
if (1.47k
partial->ma && 1.47k
partial->next->ma107
)
{89
501
89
    same = equal_on_domain(partial->ma, partial->next->ma,
502
89
          partial->dom);
503
89
    if (same < 0)
504
0
      return isl_stat_error;
505
89
    
if (89
same89
)
506
39
      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
1.47k
515
1.39k
  sol_pop_one(sol);
516
1.39k
  sol_pop_one(sol);
517
1.39k
518
1.39k
  return isl_stat_ok;
519
1.47k
}
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
5.53k
{
532
5.53k
  struct isl_partial_sol *partial;
533
5.53k
534
5.53k
  if (sol->error)
535
0
    return;
536
5.53k
537
5.53k
  partial = sol->partial;
538
5.53k
  if (!partial)
539
1.30k
    return;
540
5.53k
541
4.22k
  
if (4.22k
partial->level == 0 && 4.22k
sol->level == 01.58k
)
{1.58k
542
3.17k
    for (partial = sol->partial; 
partial3.17k
;
partial = sol->partial1.58k
)
543
1.58k
      sol_pop_one(sol);
544
1.58k
    return;
545
1.58k
  }
546
4.22k
547
2.64k
  
if (2.64k
partial->level <= sol->level2.64k
)
548
8
    return;
549
2.64k
550
2.63k
  
if (2.63k
partial->next && 2.63k
partial->next->level == partial->level1.63k
)
{1.48k
551
1.48k
    if (combine_initial_if_equal(sol) < 0)
552
0
      goto error;
553
1.48k
  } else
554
1.15k
    sol_pop_one(sol);
555
2.63k
556
2.63k
  
if (2.63k
sol->level == 02.63k
)
{999
557
1.18k
    for (partial = sol->partial; 
partial1.18k
;
partial = sol->partial181
)
558
181
      sol_pop_one(sol);
559
999
    return;
560
999
  }
561
2.63k
562
1.63k
  
if (1.63k
01.63k
)
563
0
error:    sol->error = 1;
564
1.63k
}
565
566
static void sol_dec_level(struct isl_sol *sol)
567
1.72k
{
568
1.72k
  if (sol->error)
569
0
    return;
570
1.72k
571
1.72k
  sol->level--;
572
1.72k
573
1.72k
  sol_pop(sol);
574
1.72k
}
575
576
static isl_stat sol_dec_level_wrap(struct isl_tab_callback *cb)
577
1.72k
{
578
1.72k
  struct isl_sol_callback *callback = (struct isl_sol_callback *)cb;
579
1.72k
580
1.72k
  sol_dec_level(callback->sol);
581
1.72k
582
1.72k
  return callback->sol->error ? 
isl_stat_error0
:
isl_stat_ok1.72k
;
583
1.72k
}
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
15.2k
{
593
15.2k
  struct isl_tab *tab;
594
15.2k
595
15.2k
  if (sol->error)
596
0
    return;
597
15.2k
598
15.2k
  sol->level++;
599
15.2k
  tab = sol->context->op->peek_tab(sol->context);
600
15.2k
  if (isl_tab_push_callback(tab, &sol->dec_level.callback) < 0)
601
0
    sol->error = 1;
602
15.2k
}
603
604
static void scale_rows(struct isl_mat *mat, isl_int m, int n_row)
605
9.57k
{
606
9.57k
  int i;
607
9.57k
608
9.57k
  if (isl_int_is_one(m))
609
9.56k
    return;
610
9.57k
611
51
  
for (i = 0; 18
i < n_row51
;
++i33
)
612
33
    isl_seq_scale(mat->row[i], mat->row[i], m, mat->n_col);
613
18
}
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
19.0k
{
666
19.0k
  struct isl_basic_set *bset = NULL;
667
19.0k
  struct isl_mat *mat = NULL;
668
19.0k
  unsigned off;
669
19.0k
  int row;
670
19.0k
  isl_int m;
671
19.0k
672
19.0k
  if (
sol->error || 19.0k
!tab19.0k
)
673
0
    goto error;
674
19.0k
675
19.0k
  
if (19.0k
tab->empty && 19.0k
!sol->add_empty15.0k
)
676
596
    return;
677
18.4k
  
if (18.4k
sol->context->op->is_empty(sol->context)18.4k
)
678
12.6k
    return;
679
18.4k
680
5.79k
  bset = sol_domain(sol);
681
5.79k
682
5.79k
  if (
tab->empty5.79k
)
{1.82k
683
1.82k
    sol_push_sol(sol, bset, NULL);
684
1.82k
    return;
685
1.82k
  }
686
5.79k
687
3.96k
  off = 2 + tab->M;
688
3.96k
689
3.96k
  mat = isl_mat_alloc(tab->mat->ctx, 1 + sol->n_out,
690
3.96k
              1 + tab->n_param + tab->n_div);
691
3.96k
  if (!mat)
692
0
    goto error;
693
3.96k
694
3.96k
  
isl_int_init3.96k
(m);3.96k
695
3.96k
696
3.96k
  isl_seq_clr(mat->row[0] + 1, mat->n_col - 1);
697
3.96k
  isl_int_set_si(mat->row[0][0], 1);
698
13.5k
  for (row = 0; 
row < sol->n_out13.5k
;
++row9.57k
)
{9.57k
699
9.57k
    int i = tab->n_param + row;
700
9.57k
    int r, j;
701
9.57k
702
9.57k
    isl_seq_clr(mat->row[1 + row], mat->n_col);
703
9.57k
    if (
!tab->var[i].is_row9.57k
)
{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
9.57k
710
9.57k
    r = tab->var[i].index;
711
9.57k
    if (tab->M &&
712
9.57k
        isl_int_ne(tab->mat->row[r][2], tab->mat->row[r][0]))
713
0
      isl_die(mat->ctx, isl_error_invalid,
714
9.57k
        "unbounded optimum", goto error2);
715
9.57k
    
isl_int_gcd9.57k
(m, mat->row[0][0], tab->mat->row[r][0]);9.57k
716
9.57k
    isl_int_divexact(m, tab->mat->row[r][0], m);
717
9.57k
    scale_rows(mat, m, 1 + row);
718
9.57k
    isl_int_divexact(m, mat->row[0][0], tab->mat->row[r][0]);
719
9.57k
    isl_int_mul(mat->row[1 + row][0], m, tab->mat->row[r][1]);
720
55.5k
    for (j = 0; 
j < tab->n_param55.5k
;
++j45.9k
)
{45.9k
721
45.9k
      int col;
722
45.9k
      if (tab->var[j].is_row)
723
26.9k
        continue;
724
18.9k
      col = tab->var[j].index;
725
18.9k
      isl_int_mul(mat->row[1 + row][1 + j], m,
726
18.9k
            tab->mat->row[r][off + col]);
727
18.9k
    }
728
10.8k
    for (j = 0; 
j < tab->n_div10.8k
;
++j1.27k
)
{1.27k
729
1.27k
      int col;
730
1.27k
      if (tab->var[tab->n_var - tab->n_div+j].is_row)
731
268
        continue;
732
1.01k
      col = tab->var[tab->n_var - tab->n_div+j].index;
733
1.01k
      isl_int_mul(mat->row[1 + row][1 + tab->n_param + j], m,
734
1.01k
            tab->mat->row[r][off + col]);
735
1.01k
    }
736
9.57k
    if (sol->max)
737
7.78k
      isl_seq_neg(mat->row[1 + row], mat->row[1 + row],
738
7.78k
            mat->n_col);
739
9.57k
  }
740
3.96k
741
3.96k
  
isl_int_clear3.96k
(m);3.96k
742
3.96k
743
3.96k
  sol_push_sol_mat(sol, bset, mat);
744
3.96k
  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
2.20k
{
761
2.20k
  struct isl_sol_map *sol_map = (struct isl_sol_map *) sol;
762
2.20k
  isl_map_free(sol_map->map);
763
2.20k
  isl_set_free(sol_map->empty);
764
2.20k
}
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
1.54k
{
773
1.54k
  if (
!bset || 1.54k
!sol->empty1.54k
)
774
0
    goto error;
775
1.54k
776
1.54k
  sol->empty = isl_set_grow(sol->empty, 1);
777
1.54k
  bset = isl_basic_set_simplify(bset);
778
1.54k
  bset = isl_basic_set_finalize(bset);
779
1.54k
  sol->empty = isl_set_add_basic_set(sol->empty, isl_basic_set_copy(bset));
780
1.54k
  if (!sol->empty)
781
0
    goto error;
782
1.54k
  isl_basic_set_free(bset);
783
1.54k
  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
1.54k
{
792
1.54k
  sol_map_add_empty((struct isl_sol_map *)sol, bset);
793
1.54k
}
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
2.04k
{
802
2.04k
  isl_basic_map *bmap;
803
2.04k
804
2.04k
  if (
sol->sol.error || 2.04k
!dom2.04k
||
!ma2.04k
)
805
0
    goto error;
806
2.04k
807
2.04k
  bmap = isl_basic_map_from_multi_aff2(ma, sol->sol.rational);
808
2.04k
  bmap = isl_basic_map_intersect_domain(bmap, dom);
809
2.04k
  sol->map = isl_map_grow(sol->map, 1);
810
2.04k
  sol->map = isl_map_add_basic_map(sol->map, bmap);
811
2.04k
  if (!sol->map)
812
0
    sol->sol.error = 1;
813
2.04k
  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
2.04k
{
823
2.04k
  sol_map_add((struct isl_sol_map *)sol, dom, ma);
824
2.04k
}
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
18.7k
{
836
18.7k
  int i;
837
18.7k
  unsigned off = 2 + tab->M;
838
18.7k
839
18.7k
  isl_int_set(line[0], tab->mat->row[row][1]);
840
112k
  for (i = 0; 
i < tab->n_param112k
;
++i94.1k
)
{94.1k
841
94.1k
    if (tab->var[i].is_row)
842
49.8k
      isl_int_set_si(line[1 + i], 0);
843
44.3k
    else {
844
44.3k
      int col = tab->var[i].index;
845
44.3k
      isl_int_set(line[1 + i], tab->mat->row[row][off + col]);
846
44.3k
    }
847
94.1k
  }
848
27.6k
  for (i = 0; 
i < tab->n_div27.6k
;
++i8.89k
)
{8.89k
849
8.89k
    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
5.62k
    else {
852
5.62k
      int col = tab->var[tab->n_var - tab->n_div + i].index;
853
5.62k
      isl_int_set(line[1 + tab->n_param + i],
854
5.62k
            tab->mat->row[row][off + col]);
855
5.62k
    }
856
8.89k
  }
857
18.7k
}
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
51.3k
{
867
51.3k
  int i;
868
51.3k
  unsigned off = 2 + tab->M;
869
51.3k
870
51.3k
  if (isl_int_ne(tab->mat->row[row1][1], tab->mat->row[row2][1]))
871
47.2k
    return 0;
872
51.3k
873
4.07k
  
if (4.07k
tab->M && 4.07k
isl_int_ne4.07k
(tab->mat->row[row1][2],
874
4.07k
         tab->mat->row[row2][2]))
875
1.22k
    return 0;
876
4.07k
877
6.24k
  
for (i = 0; 2.84k
i < tab->n_param + tab->n_div6.24k
;
++i3.39k
)
{6.19k
878
6.04k
    int pos = i < tab->n_param ? i :
879
153
      tab->n_var - tab->n_div + i - tab->n_param;
880
6.19k
    int col;
881
6.19k
882
6.19k
    if (tab->var[pos].is_row)
883
1.75k
      continue;
884
4.44k
    col = tab->var[pos].index;
885
4.44k
    if (isl_int_ne(tab->mat->row[row1][off + col],
886
4.44k
             tab->mat->row[row2][off + col]))
887
2.79k
      return 0;
888
4.44k
  }
889
47
  return 1;
890
2.84k
}
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
11.6k
{
899
11.6k
  struct isl_vec *ineq;
900
11.6k
901
11.6k
  ineq = isl_vec_alloc(tab->mat->ctx, 1 + tab->n_param + tab->n_div);
902
11.6k
  if (!ineq)
903
0
    return NULL;
904
11.6k
905
11.6k
  get_row_parameter_line(tab, row, ineq->el);
906
11.6k
  if (ineq)
907
11.6k
    ineq = isl_vec_normalize(ineq);
908
11.6k
909
11.6k
  return ineq;
910
11.6k
}
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
448
{
922
448
  isl_ctx *ctx = isl_vec_get_ctx(div);
923
448
  int len = div->size - 2;
924
448
925
448
  isl_seq_gcd(div->el + 2, len, &ctx->normalize_gcd);
926
448
  isl_int_gcd(ctx->normalize_gcd, ctx->normalize_gcd, div->el[0]);
927
448
928
448
  if (isl_int_is_one(ctx->normalize_gcd))
929
430
    return;
930
448
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
372
{
949
372
  struct isl_vec *div;
950
372
951
372
  div = isl_vec_alloc(tab->mat->ctx, 1 + 1 + tab->n_param + tab->n_div);
952
372
  if (!div)
953
0
    return NULL;
954
372
955
372
  
isl_int_set372
(div->el[0], tab->mat->row[row][0]);372
956
372
  get_row_parameter_line(tab, row, div->el + 1);
957
372
  isl_seq_neg(div->el + 1, div->el + 1, div->size - 1);
958
372
  normalize_div(div);
959
372
  isl_seq_fdiv_r(div->el + 1, div->el + 1, div->el[0], div->size - 1);
960
372
961
372
  return div;
962
372
}
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
76
{
977
76
  struct isl_vec *div;
978
76
979
76
  div = isl_vec_alloc(tab->mat->ctx, 1 + 1 + tab->n_param + tab->n_div);
980
76
  if (!div)
981
0
    return NULL;
982
76
983
76
  
isl_int_set76
(div->el[0], tab->mat->row[row][0]);76
984
76
  get_row_parameter_line(tab, row, div->el + 1);
985
76
  normalize_div(div);
986
76
  isl_seq_fdiv_r(div->el + 1, div->el + 1, div->el[0], div->size - 1);
987
76
988
76
  return div;
989
76
}
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 __isl_give isl_vec *ineq_for_div(__isl_keep isl_basic_set *bset,
1002
  unsigned div)
1003
76
{
1004
76
  unsigned total;
1005
76
  unsigned div_pos;
1006
76
  struct isl_vec *ineq;
1007
76
1008
76
  if (!bset)
1009
0
    return NULL;
1010
76
1011
76
  total = isl_basic_set_total_dim(bset);
1012
76
  div_pos = 1 + total - bset->n_div + div;
1013
76
1014
76
  ineq = isl_vec_alloc(bset->ctx, 1 + total);
1015
76
  if (!ineq)
1016
0
    return NULL;
1017
76
1018
76
  isl_seq_cpy(ineq->el, bset->div[div] + 1, 1 + total);
1019
76
  isl_int_neg(ineq->el[div_pos], bset->div[div][0]);
1020
76
  return ineq;
1021
76
}
1022
1023
/* Given a row in the tableau and a div that was created
1024
 * using get_row_split_div and that has been constrained to equality, i.e.,
1025
 *
1026
 *    d = floor(\sum_i {a_i} y_i) = \sum_i {a_i} y_i
1027
 *
1028
 * replace the expression "\sum_i {a_i} y_i" in the row by d,
1029
 * i.e., we subtract "\sum_i {a_i} y_i" and add 1 d.
1030
 * The coefficients of the non-parameters in the tableau have been
1031
 * verified to be integral.  We can therefore simply replace coefficient b
1032
 * by floor(b).  For the coefficients of the parameters we have
1033
 * floor(a_i) = a_i - {a_i}, while for the other coefficients, we have
1034
 * floor(b) = b.
1035
 */
1036
static struct isl_tab *set_row_cst_to_div(struct isl_tab *tab, int row, int div)
1037
76
{
1038
76
  isl_seq_fdiv_q(tab->mat->row[row] + 1, tab->mat->row[row] + 1,
1039
76
      tab->mat->row[row][0], 1 + tab->M + tab->n_col);
1040
76
1041
76
  isl_int_set_si(tab->mat->row[row][0], 1);
1042
76
1043
76
  if (
tab->var[tab->n_var - tab->n_div + div].is_row76
)
{0
1044
0
    int drow = tab->var[tab->n_var - tab->n_div + div].index;
1045
0
1046
0
    isl_assert(tab->mat->ctx,
1047
0
      isl_int_is_one(tab->mat->row[drow][0]), goto error);
1048
0
    isl_seq_combine(tab->mat->row[row] + 1,
1049
0
      tab->mat->ctx->one, tab->mat->row[row] + 1,
1050
0
      tab->mat->ctx->one, tab->mat->row[drow] + 1,
1051
0
      1 + tab->M + tab->n_col);
1052
76
  } else {
1053
76
    int dcol = tab->var[tab->n_var - tab->n_div + div].index;
1054
76
1055
76
    isl_int_add_ui(tab->mat->row[row][2 + tab->M + dcol],
1056
76
        tab->mat->row[row][2 + tab->M + dcol], 1);
1057
76
  }
1058
76
1059
76
  return tab;
1060
0
error:
1061
0
  isl_tab_free(tab);
1062
0
  return NULL;
1063
76
}
1064
1065
/* Check if the (parametric) constant of the given row is obviously
1066
 * negative, meaning that we don't need to consult the context tableau.
1067
 * If there is a big parameter and its coefficient is non-zero,
1068
 * then this coefficient determines the outcome.
1069
 * Otherwise, we check whether the constant is negative and
1070
 * all non-zero coefficients of parameters are negative and
1071
 * belong to non-negative parameters.
1072
 */
1073
static int is_obviously_neg(struct isl_tab *tab, int row)
1074
169k
{
1075
169k
  int i;
1076
169k
  int col;
1077
169k
  unsigned off = 2 + tab->M;
1078
169k
1079
169k
  if (
tab->M169k
)
{59.6k
1080
59.6k
    if (isl_int_is_pos(tab->mat->row[row][2]))
1081
31.0k
      return 0;
1082
28.5k
    
if (28.5k
isl_int_is_neg28.5k
(tab->mat->row[row][2]))
1083
0
      return 1;
1084
28.5k
  }
1085
169k
1086
138k
  
if (138k
isl_int_is_nonneg138k
(tab->mat->row[row][1]))
1087
127k
    return 0;
1088
21.1k
  
for (i = 0; 11.0k
i < tab->n_param21.1k
;
++i10.0k
)
{17.2k
1089
17.2k
    /* Eliminated parameter */
1090
17.2k
    if (tab->var[i].is_row)
1091
5.49k
      continue;
1092
11.7k
    col = tab->var[i].index;
1093
11.7k
    if (isl_int_is_zero(tab->mat->row[row][off + col]))
1094
4.47k
      continue;
1095
7.23k
    
if (7.23k
!tab->var[i].is_nonneg7.23k
)
1096
6.17k
      return 0;
1097
1.05k
    
if (1.05k
isl_int_is_pos1.05k
(tab->mat->row[row][off + col]))
1098
945
      return 0;
1099
1.05k
  }
1100
4.36k
  
for (i = 0; 3.92k
i < tab->n_div4.36k
;
++i432
)
{510
1101
510
    if (tab->var[tab->n_var - tab->n_div + i].is_row)
1102
403
      continue;
1103
107
    col = tab->var[tab->n_var - tab->n_div + i].index;
1104
107
    if (isl_int_is_zero(tab->mat->row[row][off + col]))
1105
29
      continue;
1106
78
    
if (78
!tab->var[tab->n_var - tab->n_div + i].is_nonneg78
)
1107
60
      return 0;
1108
18
    
if (18
isl_int_is_pos18
(tab->mat->row[row][off + col]))
1109
18
      return 0;
1110
18
  }
1111
3.85k
  return 1;
1112
3.92k
}
1113
1114
/* Check if the (parametric) constant of the given row is obviously
1115
 * non-negative, meaning that we don't need to consult the context tableau.
1116
 * If there is a big parameter and its coefficient is non-zero,
1117
 * then this coefficient determines the outcome.
1118
 * Otherwise, we check whether the constant is non-negative and
1119
 * all non-zero coefficients of parameters are positive and
1120
 * belong to non-negative parameters.
1121
 */
1122
static int is_obviously_nonneg(struct isl_tab *tab, int row)
1123
20.3k
{
1124
20.3k
  int i;
1125
20.3k
  int col;
1126
20.3k
  unsigned off = 2 + tab->M;
1127
20.3k
1128
20.3k
  if (
tab->M20.3k
)
{20.3k
1129
20.3k
    if (isl_int_is_pos(tab->mat->row[row][2]))
1130
8.34k
      return 1;
1131
11.9k
    
if (11.9k
isl_int_is_neg11.9k
(tab->mat->row[row][2]))
1132
0
      return 0;
1133
11.9k
  }
1134
20.3k
1135
11.9k
  
if (11.9k
isl_int_is_neg11.9k
(tab->mat->row[row][1]))
1136
4.25k
    return 0;
1137
30.3k
  
for (i = 0; 7.74k
i < tab->n_param30.3k
;
++i22.6k
)
{26.1k
1138
26.1k
    /* Eliminated parameter */
1139
26.1k
    if (tab->var[i].is_row)
1140
12.1k
      continue;
1141
14.0k
    col = tab->var[i].index;
1142
14.0k
    if (isl_int_is_zero(tab->mat->row[row][off + col]))
1143
8.88k
      continue;
1144
5.16k
    
if (5.16k
!tab->var[i].is_nonneg5.16k
)
1145
947
      return 0;
1146
4.22k
    
if (4.22k
isl_int_is_neg4.22k
(tab->mat->row[row][off + col]))
1147
2.59k
      return 0;
1148
4.22k
  }
1149
8.86k
  
for (i = 0; 4.20k
i < tab->n_div8.86k
;
++i4.66k
)
{4.73k
1150
4.73k
    if (tab->var[tab->n_var - tab->n_div + i].is_row)
1151
3.89k
      continue;
1152
845
    col = tab->var[tab->n_var - tab->n_div + i].index;
1153
845
    if (isl_int_is_zero(tab->mat->row[row][off + col]))
1154
765
      continue;
1155
80
    
if (80
!tab->var[tab->n_var - tab->n_div + i].is_nonneg80
)
1156
70
      return 0;
1157
10
    
if (10
isl_int_is_neg10
(tab->mat->row[row][off + col]))
1158
5
      return 0;
1159
10
  }
1160
4.13k
  return 1;
1161
4.20k
}
1162
1163
/* Given a row r and two columns, return the column that would
1164
 * lead to the lexicographically smallest increment in the sample
1165
 * solution when leaving the basis in favor of the row.
1166
 * Pivoting with column c will increment the sample value by a non-negative
1167
 * constant times a_{V,c}/a_{r,c}, with a_{V,c} the elements of column c
1168
 * corresponding to the non-parametric variables.
1169
 * If variable v appears in a column c_v, the a_{v,c} = 1 iff c = c_v,
1170
 * with all other entries in this virtual row equal to zero.
1171
 * If variable v appears in a row, then a_{v,c} is the element in column c
1172
 * of that row.
1173
 *
1174
 * Let v be the first variable with a_{v,c1}/a_{r,c1} != a_{v,c2}/a_{r,c2}.
1175
 * Then if a_{v,c1}/a_{r,c1} < a_{v,c2}/a_{r,c2}, i.e.,
1176
 * a_{v,c2} a_{r,c1} - a_{v,c1} a_{r,c2} > 0, c1 results in the minimal
1177
 * increment.  Otherwise, it's c2.
1178
 */
1179
static int lexmin_col_pair(struct isl_tab *tab,
1180
  int row, int col1, int col2, isl_int tmp)
1181
6.46k
{
1182
6.46k
  int i;
1183
6.46k
  isl_int *tr;
1184
6.46k
1185
6.46k
  tr = tab->mat->row[row] + 2 + tab->M;
1186
6.46k
1187
32.2k
  for (i = tab->n_param; 
i < tab->n_var - tab->n_div32.2k
;
++i25.7k
)
{32.2k
1188
32.2k
    int s1, s2;
1189
32.2k
    isl_int *r;
1190
32.2k
1191
32.2k
    if (
!tab->var[i].is_row32.2k
)
{12.0k
1192
12.0k
      if (tab->var[i].index == col1)
1193
827
        return col2;
1194
11.1k
      
if (11.1k
tab->var[i].index == col211.1k
)
1195
906
        return col1;
1196
10.2k
      continue;
1197
11.1k
    }
1198
32.2k
1199
20.1k
    
if (20.1k
tab->var[i].index == row20.1k
)
1200
173
      continue;
1201
20.1k
1202
20.0k
    r = tab->mat->row[tab->var[i].index] + 2 + tab->M;
1203
20.0k
    s1 = isl_int_sgn(r[col1]);
1204
20.0k
    s2 = isl_int_sgn(r[col2]);
1205
20.0k
    if (
s1 == 0 && 20.0k
s2 == 016.4k
)
1206
14.1k
      continue;
1207
5.86k
    
if (5.86k
s1 < s25.86k
)
1208
2.53k
      return col1;
1209
3.33k
    
if (3.33k
s2 < s13.33k
)
1210
776
      return col2;
1211
3.33k
1212
2.55k
    
isl_int_mul2.55k
(tmp, r[col2], tr[col1]);2.55k
1213
2.55k
    isl_int_submul(tmp, r[col1], tr[col2]);
1214
2.55k
    if (isl_int_is_pos(tmp))
1215
901
      return col1;
1216
1.65k
    
if (1.65k
isl_int_is_neg1.65k
(tmp))
1217
518
      return col2;
1218
1.65k
  }
1219
0
  return -1;
1220
6.46k
}
1221
1222
/* Given a row in the tableau, find and return the column that would
1223
 * result in the lexicographically smallest, but positive, increment
1224
 * in the sample point.
1225
 * If there is no such column, then return tab->n_col.
1226
 * If anything goes wrong, return -1.
1227
 */
1228
static int lexmin_pivot_col(struct isl_tab *tab, int row)
1229
9.14k
{
1230
9.14k
  int j;
1231
9.14k
  int col = tab->n_col;
1232
9.14k
  isl_int *tr;
1233
9.14k
  isl_int tmp;
1234
9.14k
1235
9.14k
  tr = tab->mat->row[row] + 2 + tab->M;
1236
9.14k
1237
9.14k
  isl_int_init(tmp);
1238
9.14k
1239
79.4k
  for (j = tab->n_dead; 
j < tab->n_col79.4k
;
++j70.3k
)
{70.3k
1240
70.3k
    if (tab->col_var[j] >= 0 &&
1241
51.2k
        (tab->col_var[j] < tab->n_param  ||
1242
40.2k
        tab->col_var[j] >= tab->n_var - tab->n_div))
1243
13.1k
      continue;
1244
70.3k
1245
57.2k
    
if (57.2k
!57.2k
isl_int_is_pos57.2k
(tr[j]))
1246
43.8k
      continue;
1247
57.2k
1248
13.3k
    
if (13.3k
col == tab->n_col13.3k
)
1249
6.91k
      col = j;
1250
13.3k
    else
1251
6.46k
      col = lexmin_col_pair(tab, row, col, j, tmp);
1252
13.3k
    isl_assert(tab->mat->ctx, col >= 0, goto error);
1253
13.3k
  }
1254
9.14k
1255
9.14k
  
isl_int_clear9.14k
(tmp);9.14k
1256
9.14k
  return col;
1257
0
error:
1258
0
  isl_int_clear(tmp);
1259
0
  return -1;
1260
9.14k
}
1261
1262
/* Return the first known violated constraint, i.e., a non-negative
1263
 * constraint that currently has an either obviously negative value
1264
 * or a previously determined to be negative value.
1265
 *
1266
 * If any constraint has a negative coefficient for the big parameter,
1267
 * if any, then we return one of these first.
1268
 */
1269
static int first_neg(struct isl_tab *tab)
1270
28.1k
{
1271
28.1k
  int row;
1272
28.1k
1273
28.1k
  if (tab->M)
1274
174k
    
for (row = tab->n_redundant; 19.9k
row < tab->n_row174k
;
++row155k
)
{157k
1275
157k
      if (!isl_tab_var_from_row(tab, row)->is_nonneg)
1276
40.8k
        continue;
1277
116k
      
if (116k
!116k
isl_int_is_neg116k
(tab->mat->row[row][2]))
1278
114k
        continue;
1279
2.53k
      
if (2.53k
tab->row_sign2.53k
)
1280
2.53k
        tab->row_sign[row] = isl_tab_row_neg;
1281
2.53k
      return row;
1282
116k
    }
1283
267k
  
for (row = tab->n_redundant; 25.6k
row < tab->n_row267k
;
++row241k
)
{248k
1284
248k
    if (!isl_tab_var_from_row(tab, row)->is_nonneg)
1285
35.1k
      continue;
1286
213k
    
if (213k
tab->row_sign213k
)
{103k
1287
103k
      if (tab->row_sign[row] == 0 &&
1288
59.6k
          is_obviously_neg(tab, row))
1289
203
        tab->row_sign[row] = isl_tab_row_neg;
1290
103k
      if (tab->row_sign[row] != isl_tab_row_neg)
1291
100k
        continue;
1292
110k
    } else 
if (110k
!is_obviously_neg(tab, row)110k
)
1293
106k
      continue;
1294
6.61k
    return row;
1295
213k
  }
1296
19.0k
  return -1;
1297
25.6k
}
1298
1299
/* Check whether the invariant that all columns are lexico-positive
1300
 * is satisfied.  This function is not called from the current code
1301
 * but is useful during debugging.
1302
 */
1303
static void check_lexpos(struct isl_tab *tab) __attribute__ ((unused));
1304
static void check_lexpos(struct isl_tab *tab)
1305
0
{
1306
0
  unsigned off = 2 + tab->M;
1307
0
  int col;
1308
0
  int var;
1309
0
  int row;
1310
0
1311
0
  for (col = tab->n_dead; col < tab->n_col; ++col) {
1312
0
    if (tab->col_var[col] >= 0 &&
1313
0
        (tab->col_var[col] < tab->n_param ||
1314
0
         tab->col_var[col] >= tab->n_var - tab->n_div))
1315
0
      continue;
1316
0
    for (var = tab->n_param; var < tab->n_var - tab->n_div; ++var) {
1317
0
      if (!tab->var[var].is_row) {
1318
0
        if (tab->var[var].index == col)
1319
0
          break;
1320
0
        else
1321
0
          continue;
1322
0
      }
1323
0
      row = tab->var[var].index;
1324
0
      if (isl_int_is_zero(tab->mat->row[row][off + col]))
1325
0
        continue;
1326
0
      if (isl_int_is_pos(tab->mat->row[row][off + col]))
1327
0
        break;
1328
0
      fprintf(stderr, "lexneg column %d (row %d)\n",
1329
0
        col, row);
1330
0
    }
1331
0
    if (var >= tab->n_var - tab->n_div)
1332
0
      fprintf(stderr, "zero column %d\n", col);
1333
0
  }
1334
0
}
1335
1336
/* Report to the caller that the given constraint is part of an encountered
1337
 * conflict.
1338
 */
1339
static int report_conflicting_constraint(struct isl_tab *tab, int con)
1340
1.17k
{
1341
1.17k
  return tab->conflict(con, tab->conflict_user);
1342
1.17k
}
1343
1344
/* Given a conflicting row in the tableau, report all constraints
1345
 * involved in the row to the caller.  That is, the row itself
1346
 * (if it represents a constraint) and all constraint columns with
1347
 * non-zero (and therefore negative) coefficients.
1348
 */
1349
static int report_conflict(struct isl_tab *tab, int row)
1350
2.22k
{
1351
2.22k
  int j;
1352
2.22k
  isl_int *tr;
1353
2.22k
1354
2.22k
  if (!tab->conflict)
1355
1.79k
    return 0;
1356
2.22k
1357
436
  
if (436
tab->row_var[row] < 0 &&436
1358
436
      report_conflicting_constraint(tab, ~tab->row_var[row]) < 0)
1359
0
    return -1;
1360
436
1361
436
  tr = tab->mat->row[row] + 2 + tab->M;
1362
436
1363
5.17k
  for (j = tab->n_dead; 
j < tab->n_col5.17k
;
++j4.73k
)
{4.73k
1364
4.73k
    if (tab->col_var[j] >= 0 &&
1365
3.80k
        (tab->col_var[j] < tab->n_param  ||
1366
3.80k
        tab->col_var[j] >= tab->n_var - tab->n_div))
1367
0
      continue;
1368
4.73k
1369
4.73k
    
if (4.73k
!4.73k
isl_int_is_neg4.73k
(tr[j]))
1370
3.93k
      continue;
1371
4.73k
1372
798
    
if (798
tab->col_var[j] < 0 &&798
1373
737
        report_conflicting_constraint(tab, ~tab->col_var[j]) < 0)
1374
0
      return -1;
1375
798
  }
1376
436
1377
436
  return 0;
1378
436
}
1379
1380
/* Resolve all known or obviously violated constraints through pivoting.
1381
 * In particular, as long as we can find any violated constraint, we
1382
 * look for a pivoting column that would result in the lexicographically
1383
 * smallest increment in the sample point.  If there is no such column
1384
 * then the tableau is infeasible.
1385
 */
1386
static int restore_lexmin(struct isl_tab *tab) WARN_UNUSED;
1387
static int restore_lexmin(struct isl_tab *tab)
1388
21.2k
{
1389
21.2k
  int row, col;
1390
21.2k
1391
21.2k
  if (!tab)
1392
0
    return -1;
1393
21.2k
  
if (21.2k
tab->empty21.2k
)
1394
32
    return 0;
1395
28.1k
  
while (21.2k
(row = first_neg(tab)) != -128.1k
)
{9.14k
1396
9.14k
    col = lexmin_pivot_col(tab, row);
1397
9.14k
    if (
col >= tab->n_col9.14k
)
{2.22k
1398
2.22k
      if (report_conflict(tab, row) < 0)
1399
0
        return -1;
1400
2.22k
      
if (2.22k
isl_tab_mark_empty(tab) < 02.22k
)
1401
0
        return -1;
1402
2.22k
      return 0;
1403
2.22k
    }
1404
6.91k
    
if (6.91k
col < 06.91k
)
1405
0
      return -1;
1406
6.91k
    
if (6.91k
isl_tab_pivot(tab, row, col) < 06.91k
)
1407
0
      return -1;
1408
6.91k
  }
1409
19.0k
  return 0;
1410
21.2k
}
1411
1412
/* Given a row that represents an equality, look for an appropriate
1413
 * pivoting column.
1414
 * In particular, if there are any non-zero coefficients among
1415
 * the non-parameter variables, then we take the last of these
1416
 * variables.  Eliminating this variable in terms of the other
1417
 * variables and/or parameters does not influence the property
1418
 * that all column in the initial tableau are lexicographically
1419
 * positive.  The row corresponding to the eliminated variable
1420
 * will only have non-zero entries below the diagonal of the
1421
 * initial tableau.  That is, we transform
1422
 *
1423
 *    I       I
1424
 *      1   into    a
1425
 *        I         I
1426
 *
1427
 * If there is no such non-parameter variable, then we are dealing with
1428
 * pure parameter equality and we pick any parameter with coefficient 1 or -1
1429
 * for elimination.  This will ensure that the eliminated parameter
1430
 * always has an integer value whenever all the other parameters are integral.
1431
 * If there is no such parameter then we return -1.
1432
 */
1433
static int last_var_col_or_int_par_col(struct isl_tab *tab, int row)
1434
15.1k
{
1435
15.1k
  unsigned off = 2 + tab->M;
1436
15.1k
  int i;
1437
15.1k
1438
59.8k
  for (i = tab->n_var - tab->n_div - 1; 
i >= 0 && 59.8k
i >= tab->n_param59.2k
;
--i44.7k
)
{52.6k
1439
52.6k
    int col;
1440
52.6k
    if (tab->var[i].is_row)
1441
30.3k
      continue;
1442
22.3k
    col = tab->var[i].index;
1443
22.3k
    if (col <= tab->n_dead)
1444
1.68k
      continue;
1445
20.6k
    
if (20.6k
!20.6k
isl_int_is_zero20.6k
(tab->mat->row[row][off + col]))
1446
7.90k
      return col;
1447
20.6k
  }
1448
21.9k
  
for (i = tab->n_dead; 7.21k
i < tab->n_col21.9k
;
++i14.6k
)
{21.7k
1449
21.7k
    if (isl_int_is_one(tab->mat->row[row][off + i]))
1450
3.19k
      return i;
1451
18.5k
    
if (18.5k
isl_int_is_negone18.5k
(tab->mat->row[row][off + i]))
1452
3.85k
      return i;
1453
18.5k
  }
1454
169
  return -1;
1455
7.21k
}
1456
1457
/* Add an equality that is known to be valid to the tableau.
1458
 * We first check if we can eliminate a variable or a parameter.
1459
 * If not, we add the equality as two inequalities.
1460
 * In this case, the equality was a pure parameter equality and there
1461
 * is no need to resolve any constraint violations.
1462
 *
1463
 * This function assumes that at least two more rows and at least
1464
 * two more elements in the constraint array are available in the tableau.
1465
 */
1466
static struct isl_tab *add_lexmin_valid_eq(struct isl_tab *tab, isl_int *eq)
1467
15.1k
{
1468
15.1k
  int i;
1469
15.1k
  int r;
1470
15.1k
1471
15.1k
  if (!tab)
1472
0
    return NULL;
1473
15.1k
  r = isl_tab_add_row(tab, eq);
1474
15.1k
  if (r < 0)
1475
0
    goto error;
1476
15.1k
1477
15.1k
  r = tab->con[r].index;
1478
15.1k
  i = last_var_col_or_int_par_col(tab, r);
1479
15.1k
  if (
i < 015.1k
)
{169
1480
169
    tab->con[r].is_nonneg = 1;
1481
169
    if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r]) < 0)
1482
0
      goto error;
1483
169
    isl_seq_neg(eq, eq, 1 + tab->n_var);
1484
169
    r = isl_tab_add_row(tab, eq);
1485
169
    if (r < 0)
1486
0
      goto error;
1487
169
    tab->con[r].is_nonneg = 1;
1488
169
    if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r]) < 0)
1489
0
      goto error;
1490
14.9k
  } else {
1491
14.9k
    if (isl_tab_pivot(tab, r, i) < 0)
1492
0
      goto error;
1493
14.9k
    
if (14.9k
isl_tab_kill_col(tab, i) < 014.9k
)
1494
0
      goto error;
1495
14.9k
    tab->n_eq++;
1496
14.9k
  }
1497
15.1k
1498
15.1k
  return tab;
1499
0
error:
1500
0
  isl_tab_free(tab);
1501
0
  return NULL;
1502
15.1k
}
1503
1504
/* Check if the given row is a pure constant.
1505
 */
1506
static int is_constant(struct isl_tab *tab, int row)
1507
353
{
1508
353
  unsigned off = 2 + tab->M;
1509
353
1510
353
  return isl_seq_first_non_zero(tab->mat->row[row] + off + tab->n_dead,
1511
353
          tab->n_col - tab->n_dead) == -1;
1512
353
}
1513
1514
/* Add an equality that may or may not be valid to the tableau.
1515
 * If the resulting row is a pure constant, then it must be zero.
1516
 * Otherwise, the resulting tableau is empty.
1517
 *
1518
 * If the row is not a pure constant, then we add two inequalities,
1519
 * each time checking that they can be satisfied.
1520
 * In the end we try to use one of the two constraints to eliminate
1521
 * a column.
1522
 *
1523
 * This function assumes that at least two more rows and at least
1524
 * two more elements in the constraint array are available in the tableau.
1525
 */
1526
static int add_lexmin_eq(struct isl_tab *tab, isl_int *eq) WARN_UNUSED;
1527
static int add_lexmin_eq(struct isl_tab *tab, isl_int *eq)
1528
353
{
1529
353
  int r1, r2;
1530
353
  int row;
1531
353
  struct isl_tab_undo *snap;
1532
353
1533
353
  if (!tab)
1534
0
    return -1;
1535
353
  snap = isl_tab_snap(tab);
1536
353
  r1 = isl_tab_add_row(tab, eq);
1537
353
  if (r1 < 0)
1538
0
    return -1;
1539
353
  tab->con[r1].is_nonneg = 1;
1540
353
  if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r1]) < 0)
1541
0
    return -1;
1542
353
1543
353
  row = tab->con[r1].index;
1544
353
  if (
is_constant(tab, row)353
)
{51
1545
51
    if (
!51
isl_int_is_zero51
(tab->mat->row[row][1]) ||
1546
51
        
(tab->M && 51
!0
isl_int_is_zero0
(tab->mat->row[row][2])))
{0
1547
0
      if (isl_tab_mark_empty(tab) < 0)
1548
0
        return -1;
1549
0
      return 0;
1550
0
    }
1551
51
    
if (51
isl_tab_rollback(tab, snap) < 051
)
1552
0
      return -1;
1553
51
    return 0;
1554
51
  }
1555
353
1556
302
  
if (302
restore_lexmin(tab) < 0302
)
1557
0
    return -1;
1558
302
  
if (302
tab->empty302
)
1559
20
    return 0;
1560
302
1561
282
  isl_seq_neg(eq, eq, 1 + tab->n_var);
1562
282
1563
282
  r2 = isl_tab_add_row(tab, eq);
1564
282
  if (r2 < 0)
1565
0
    return -1;
1566
282
  tab->con[r2].is_nonneg = 1;
1567
282
  if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r2]) < 0)
1568
0
    return -1;
1569
282
1570
282
  
if (282
restore_lexmin(tab) < 0282
)
1571
0
    return -1;
1572
282
  
if (282
tab->empty282
)
1573
0
    return 0;
1574
282
1575
282
  
if (282
!tab->con[r1].is_row282
)
{0
1576
0
    if (isl_tab_kill_col(tab, tab->con[r1].index) < 0)
1577
0
      return -1;
1578
282
  } else 
if (282
!tab->con[r2].is_row282
)
{0
1579
0
    if (isl_tab_kill_col(tab, tab->con[r2].index) < 0)
1580
0
      return -1;
1581
0
  }
1582
282
1583
282
  
if (282
tab->bmap282
)
{0
1584
0
    tab->bmap = isl_basic_map_add_ineq(tab->bmap, eq);
1585
0
    if (isl_tab_push(tab, isl_tab_undo_bmap_ineq) < 0)
1586
0
      return -1;
1587
0
    isl_seq_neg(eq, eq, 1 + tab->n_var);
1588
0
    tab->bmap = isl_basic_map_add_ineq(tab->bmap, eq);
1589
0
    isl_seq_neg(eq, eq, 1 + tab->n_var);
1590
0
    if (isl_tab_push(tab, isl_tab_undo_bmap_ineq) < 0)
1591
0
      return -1;
1592
0
    
if (0
!tab->bmap0
)
1593
0
      return -1;
1594
0
  }
1595
282
1596
282
  return 0;
1597
282
}
1598
1599
/* Add an inequality to the tableau, resolving violations using
1600
 * restore_lexmin.
1601
 *
1602
 * This function assumes that at least one more row and at least
1603
 * one more element in the constraint array are available in the tableau.
1604
 */
1605
static struct isl_tab *add_lexmin_ineq(struct isl_tab *tab, isl_int *ineq)
1606
14.4k
{
1607
14.4k
  int r;
1608
14.4k
1609
14.4k
  if (!tab)
1610
0
    return NULL;
1611
14.4k
  
if (14.4k
tab->bmap14.4k
)
{0
1612
0
    tab->bmap = isl_basic_map_add_ineq(tab->bmap, ineq);
1613
0
    if (isl_tab_push(tab, isl_tab_undo_bmap_ineq) < 0)
1614
0
      goto error;
1615
0
    
if (0
!tab->bmap0
)
1616
0
      goto error;
1617
0
  }
1618
14.4k
  r = isl_tab_add_row(tab, ineq);
1619
14.4k
  if (r < 0)
1620
0
    goto error;
1621
14.4k
  tab->con[r].is_nonneg = 1;
1622
14.4k
  if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r]) < 0)
1623
0
    goto error;
1624
14.4k
  
if (14.4k
isl_tab_row_is_redundant(tab, tab->con[r].index)14.4k
)
{302
1625
302
    if (isl_tab_mark_redundant(tab, tab->con[r].index) < 0)
1626
0
      goto error;
1627
302
    return tab;
1628
302
  }
1629
14.4k
1630
14.1k
  
if (14.1k
restore_lexmin(tab) < 014.1k
)
1631
0
    goto error;
1632
14.1k
  
if (14.1k
!tab->empty && 14.1k
tab->con[r].is_row13.7k
&&
1633
10.8k
     isl_tab_row_is_redundant(tab, tab->con[r].index))
1634
0
    
if (0
isl_tab_mark_redundant(tab, tab->con[r].index) < 00
)
1635
0
      goto error;
1636
14.1k
  return tab;
1637
0
error:
1638
0
  isl_tab_free(tab);
1639
0
  return NULL;
1640
14.1k
}
1641
1642
/* Check if the coefficients of the parameters are all integral.
1643
 */
1644
static int integer_parameter(struct isl_tab *tab, int row)
1645
14.5k
{
1646
14.5k
  int i;
1647
14.5k
  int col;
1648
14.5k
  unsigned off = 2 + tab->M;
1649
14.5k
1650
64.8k
  for (i = 0; 
i < tab->n_param64.8k
;
++i50.2k
)
{50.6k
1651
50.6k
    /* Eliminated parameter */
1652
50.6k
    if (tab->var[i].is_row)
1653
27.9k
      continue;
1654
22.7k
    col = tab->var[i].index;
1655
22.7k
    if (
!22.7k
isl_int_is_divisible_by22.7k
(tab->mat->row[row][off + col],
1656
22.7k
            tab->mat->row[row][0]))
1657
433
      return 0;
1658
22.7k
  }
1659
19.5k
  
for (i = 0; 14.1k
i < tab->n_div19.5k
;
++i5.38k
)
{5.40k
1660
5.40k
    if (tab->var[tab->n_var - tab->n_div + i].is_row)
1661
1.98k
      continue;
1662
3.41k
    col = tab->var[tab->n_var - tab->n_div + i].index;
1663
3.41k
    if (
!3.41k
isl_int_is_divisible_by3.41k
(tab->mat->row[row][off + col],
1664
3.41k
            tab->mat->row[row][0]))
1665
15
      return 0;
1666
3.41k
  }
1667
14.1k
  return 1;
1668
14.1k
}
1669
1670
/* Check if the coefficients of the non-parameter variables are all integral.
1671
 */
1672
static int integer_variable(struct isl_tab *tab, int row)
1673
861
{
1674
861
  int i;
1675
861
  unsigned off = 2 + tab->M;
1676
861
1677
2.10k
  for (i = tab->n_dead; 
i < tab->n_col2.10k
;
++i1.23k
)
{2.02k
1678
2.02k
    if (tab->col_var[i] >= 0 &&
1679
880
        (tab->col_var[i] < tab->n_param ||
1680
133
         tab->col_var[i] >= tab->n_var - tab->n_div))
1681
852
      continue;
1682
1.17k
    
if (1.17k
!1.17k
isl_int_is_divisible_by1.17k
(tab->mat->row[row][off + i],
1683
1.17k
            tab->mat->row[row][0]))
1684
785
      return 0;
1685
1.17k
  }
1686
76
  return 1;
1687
861
}
1688
1689
/* Check if the constant term is integral.
1690
 */
1691
static int integer_constant(struct isl_tab *tab, int row)
1692
14.5k
{
1693
14.5k
  return isl_int_is_divisible_by(tab->mat->row[row][1],
1694
14.5k
          tab->mat->row[row][0]);
1695
14.5k
}
1696
1697
#define I_CST 1 << 0
1698
#define I_PAR 1 << 1
1699
#define I_VAR 1 << 2
1700
1701
/* Check for next (non-parameter) variable after "var" (first if var == -1)
1702
 * that is non-integer and therefore requires a cut and return
1703
 * the index of the variable.
1704
 * For parametric tableaus, there are three parts in a row,
1705
 * the constant, the coefficients of the parameters and the rest.
1706
 * For each part, we check whether the coefficients in that part
1707
 * are all integral and if so, set the corresponding flag in *f.
1708
 * If the constant and the parameter part are integral, then the
1709
 * current sample value is integral and no cut is required
1710
 * (irrespective of whether the variable part is integral).
1711
 */
1712
static int next_non_integer_var(struct isl_tab *tab, int var, int *f)
1713
5.49k
{
1714
5.49k
  var = var < 0 ? 
tab->n_param5.49k
:
var + 10
;
1715
5.49k
1716
25.1k
  for (; 
var < tab->n_var - tab->n_div25.1k
;
++var19.6k
)
{20.5k
1717
20.5k
    int flags = 0;
1718
20.5k
    int row;
1719
20.5k
    if (!tab->var[var].is_row)
1720
5.94k
      continue;
1721
14.5k
    row = tab->var[var].index;
1722
14.5k
    if (integer_constant(tab, row))
1723
13.9k
      ISL_FL_SET(flags, I_CST);
1724
14.5k
    if (integer_parameter(tab, row))
1725
14.1k
      ISL_FL_SET(flags, I_PAR);
1726
14.5k
    if (
ISL_FL_ISSET14.5k
(flags, I_CST) && 14.5k
ISL_FL_ISSET13.9k
(flags, I_PAR))
1727
13.7k
      continue;
1728
861
    
if (861
integer_variable(tab, row)861
)
1729
76
      ISL_FL_SET(flags, I_VAR);
1730
861
    *f = flags;
1731
861
    return var;
1732
14.5k
  }
1733
4.63k
  return -1;
1734
5.49k
}
1735
1736
/* Check for first (non-parameter) variable that is non-integer and
1737
 * therefore requires a cut and return the corresponding row.
1738
 * For parametric tableaus, there are three parts in a row,
1739
 * the constant, the coefficients of the parameters and the rest.
1740
 * For each part, we check whether the coefficients in that part
1741
 * are all integral and if so, set the corresponding flag in *f.
1742
 * If the constant and the parameter part are integral, then the
1743
 * current sample value is integral and no cut is required
1744
 * (irrespective of whether the variable part is integral).
1745
 */
1746
static int first_non_integer_row(struct isl_tab *tab, int *f)
1747
4.74k
{
1748
4.74k
  int var = next_non_integer_var(tab, -1, f);
1749
4.74k
1750
3.96k
  return var < 0 ? 
-13.96k
:
tab->var[var].index783
;
1751
4.74k
}
1752
1753
/* Add a (non-parametric) cut to cut away the non-integral sample
1754
 * value of the given row.
1755
 *
1756
 * If the row is given by
1757
 *
1758
 *  m r = f + \sum_i a_i y_i
1759
 *
1760
 * then the cut is
1761
 *
1762
 *  c = - {-f/m} + \sum_i {a_i/m} y_i >= 0
1763
 *
1764
 * The big parameter, if any, is ignored, since it is assumed to be big
1765
 * enough to be divisible by any integer.
1766
 * If the tableau is actually a parametric tableau, then this function
1767
 * is only called when all coefficients of the parameters are integral.
1768
 * The cut therefore has zero coefficients for the parameters.
1769
 *
1770
 * The current value is known to be negative, so row_sign, if it
1771
 * exists, is set accordingly.
1772
 *
1773
 * Return the row of the cut or -1.
1774
 */
1775
static int add_cut(struct isl_tab *tab, int row)
1776
413
{
1777
413
  int i;
1778
413
  int r;
1779
413
  isl_int *r_row;
1780
413
  unsigned off = 2 + tab->M;
1781
413
1782
413
  if (isl_tab_extend_cons(tab, 1) < 0)
1783
0
    return -1;
1784
413
  r = isl_tab_allocate_con(tab);
1785
413
  if (r < 0)
1786
0
    return -1;
1787
413
1788
413
  r_row = tab->mat->row[tab->con[r].index];
1789
413
  isl_int_set(r_row[0], tab->mat->row[row][0]);
1790
413
  isl_int_neg(r_row[1], tab->mat->row[row][1]);
1791
413
  isl_int_fdiv_r(r_row[1], r_row[1], tab->mat->row[row][0]);
1792
413
  isl_int_neg(r_row[1], r_row[1]);
1793
413
  if (tab->M)
1794
335
    isl_int_set_si(r_row[2], 0);
1795
4.46k
  for (i = 0; 
i < tab->n_col4.46k
;
++i4.05k
)
1796
4.05k
    isl_int_fdiv_r(r_row[off + i],
1797
413
      tab->mat->row[row][off + i], tab->mat->row[row][0]);
1798
413
1799
413
  tab->con[r].is_nonneg = 1;
1800
413
  if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r]) < 0)
1801
0
    return -1;
1802
413
  
if (413
tab->row_sign413
)
1803
335
    tab->row_sign[tab->con[r].index] = isl_tab_row_neg;
1804
413
1805
413
  return tab->con[r].index;
1806
413
}
1807
1808
0
#define CUT_ALL 1
1809
1.19k
#define CUT_ONE 0
1810
1811
/* Given a non-parametric tableau, add cuts until an integer
1812
 * sample point is obtained or until the tableau is determined
1813
 * to be integer infeasible.
1814
 * As long as there is any non-integer value in the sample point,
1815
 * we add appropriate cuts, if possible, for each of these
1816
 * non-integer values and then resolve the violated
1817
 * cut constraints using restore_lexmin.
1818
 * If one of the corresponding rows is equal to an integral
1819
 * combination of variables/constraints plus a non-integral constant,
1820
 * then there is no way to obtain an integer point and we return
1821
 * a tableau that is marked empty.
1822
 * The parameter cutting_strategy controls the strategy used when adding cuts
1823
 * to remove non-integer points. CUT_ALL adds all possible cuts
1824
 * before continuing the search. CUT_ONE adds only one cut at a time.
1825
 */
1826
static struct isl_tab *cut_to_integer_lexmin(struct isl_tab *tab,
1827
  int cutting_strategy)
1828
1.11k
{
1829
1.11k
  int var;
1830
1.11k
  int row;
1831
1.11k
  int flags;
1832
1.11k
1833
1.11k
  if (!tab)
1834
0
    return NULL;
1835
1.11k
  
if (1.11k
tab->empty1.11k
)
1836
446
    return tab;
1837
1.11k
1838
748
  
while (672
(var = next_non_integer_var(tab, -1, &flags)) != -1748
)
{78
1839
78
    do {
1840
78
      if (
ISL_FL_ISSET78
(flags, I_VAR))
{0
1841
0
        if (isl_tab_mark_empty(tab) < 0)
1842
0
          goto error;
1843
0
        return tab;
1844
0
      }
1845
78
      row = tab->var[var].index;
1846
78
      row = add_cut(tab, row);
1847
78
      if (row < 0)
1848
0
        goto error;
1849
78
      
if (78
cutting_strategy == 78
CUT_ONE78
)
1850
78
        break;
1851
0
    } while ((var = next_non_integer_var(tab, var, &flags)) != -1);
1852
78
    
if (78
restore_lexmin(tab) < 078
)
1853
0
      goto error;
1854
78
    
if (78
tab->empty78
)
1855
2
      break;
1856
78
  }
1857
672
  return tab;
1858
0
error:
1859
0
  isl_tab_free(tab);
1860
0
  return NULL;
1861
672
}
1862
1863
/* Check whether all the currently active samples also satisfy the inequality
1864
 * "ineq" (treated as an equality if eq is set).
1865
 * Remove those samples that do not.
1866
 */
1867
static struct isl_tab *check_samples(struct isl_tab *tab, isl_int *ineq, int eq)
1868
10.5k
{
1869
10.5k
  int i;
1870
10.5k
  isl_int v;
1871
10.5k
1872
10.5k
  if (!tab)
1873
0
    return NULL;
1874
10.5k
1875
10.5k
  
isl_assert10.5k
(tab->mat->ctx, tab->bmap, goto error);10.5k
1876
10.5k
  
isl_assert10.5k
(tab->mat->ctx, tab->samples, goto error);10.5k
1877
10.5k
  
isl_assert10.5k
(tab->mat->ctx, tab->samples->n_col == 1 + tab->n_var, goto error);10.5k
1878
10.5k
1879
10.5k
  
isl_int_init10.5k
(v);10.5k
1880
26.1k
  for (i = tab->n_outside; 
i < tab->n_sample26.1k
;
++i15.5k
)
{15.5k
1881
15.5k
    int sgn;
1882
15.5k
    isl_seq_inner_product(ineq, tab->samples->row[i],
1883
15.5k
          1 + tab->n_var, &v);
1884
15.5k
    sgn = isl_int_sgn(v);
1885
15.5k
    if (
eq ? 15.5k
(sgn == 0)6.69k
:
(sgn >= 0)8.83k
)
1886
11.0k
      continue;
1887
4.46k
    tab = isl_tab_drop_sample(tab, i);
1888
4.46k
    if (!tab)
1889
0
      break;
1890
4.46k
  }
1891
10.5k
  isl_int_clear(v);
1892
10.5k
1893
10.5k
  return tab;
1894
0
error:
1895
0
  isl_tab_free(tab);
1896
0
  return NULL;
1897
10.5k
}
1898
1899
/* Check whether the sample value of the tableau is finite,
1900
 * i.e., either the tableau does not use a big parameter, or
1901
 * all values of the variables are equal to the big parameter plus
1902
 * some constant.  This constant is the actual sample value.
1903
 */
1904
static int sample_is_finite(struct isl_tab *tab)
1905
0
{
1906
0
  int i;
1907
0
1908
0
  if (!tab->M)
1909
0
    return 1;
1910
0
1911
0
  
for (i = 0; 0
i < tab->n_var0
;
++i0
)
{0
1912
0
    int row;
1913
0
    if (!tab->var[i].is_row)
1914
0
      return 0;
1915
0
    row = tab->var[i].index;
1916
0
    if (isl_int_ne(tab->mat->row[row][0], tab->mat->row[row][2]))
1917
0
      return 0;
1918
0
  }
1919
0
  return 1;
1920
0
}
1921
1922
/* Check if the context tableau of sol has any integer points.
1923
 * Leave tab in empty state if no integer point can be found.
1924
 * If an integer point can be found and if moreover it is finite,
1925
 * then it is added to the list of sample values.
1926
 *
1927
 * This function is only called when none of the currently active sample
1928
 * values satisfies the most recently added constraint.
1929
 */
1930
static struct isl_tab *check_integer_feasible(struct isl_tab *tab)
1931
0
{
1932
0
  struct isl_tab_undo *snap;
1933
0
1934
0
  if (!tab)
1935
0
    return NULL;
1936
0
1937
0
  snap = isl_tab_snap(tab);
1938
0
  if (isl_tab_push_basis(tab) < 0)
1939
0
    goto error;
1940
0
1941
0
  
tab = cut_to_integer_lexmin(tab, 0
CUT_ALL0
);
1942
0
  if (!tab)
1943
0
    goto error;
1944
0
1945
0
  
if (0
!tab->empty && 0
sample_is_finite(tab)0
)
{0
1946
0
    struct isl_vec *sample;
1947
0
1948
0
    sample = isl_tab_get_sample_value(tab);
1949
0
1950
0
    if (isl_tab_add_sample(tab, sample) < 0)
1951
0
      goto error;
1952
0
  }
1953
0
1954
0
  
if (0
!tab->empty && 0
isl_tab_rollback(tab, snap) < 00
)
1955
0
    goto error;
1956
0
1957
0
  return tab;
1958
0
error:
1959
0
  isl_tab_free(tab);
1960
0
  return NULL;
1961
0
}
1962
1963
/* Check if any of the currently active sample values satisfies
1964
 * the inequality "ineq" (an equality if eq is set).
1965
 */
1966
static int tab_has_valid_sample(struct isl_tab *tab, isl_int *ineq, int eq)
1967
19.9k
{
1968
19.9k
  int i;
1969
19.9k
  isl_int v;
1970
19.9k
1971
19.9k
  if (!tab)
1972
0
    return -1;
1973
19.9k
1974
19.9k
  
isl_assert19.9k
(tab->mat->ctx, tab->bmap, return -1);19.9k
1975
19.9k
  
isl_assert19.9k
(tab->mat->ctx, tab->samples, return -1);19.9k
1976
19.9k
  
isl_assert19.9k
(tab->mat->ctx, tab->samples->n_col == 1 + tab->n_var, return -1);19.9k
1977
19.9k
1978
19.9k
  
isl_int_init19.9k
(v);19.9k
1979
33.2k
  for (i = tab->n_outside; 
i < tab->n_sample33.2k
;
++i13.2k
)
{19.9k
1980
19.9k
    int sgn;
1981
19.9k
    isl_seq_inner_product(ineq, tab->samples->row[i],
1982
19.9k
          1 + tab->n_var, &v);
1983
19.9k
    sgn = isl_int_sgn(v);
1984
19.9k
    if (
eq ? 19.9k
(sgn == 0)6.60k
:
(sgn >= 0)13.3k
)
1985
6.68k
      break;
1986
19.9k
  }
1987
19.9k
  isl_int_clear(v);
1988
19.9k
1989
19.9k
  return i < tab->n_sample;
1990
19.9k
}
1991
1992
/* Insert a div specified by "div" to the tableau "tab" at position "pos" and
1993
 * return isl_bool_true if the div is obviously non-negative.
1994
 */
1995
static isl_bool context_tab_insert_div(struct isl_tab *tab, int pos,
1996
  __isl_keep isl_vec *div,
1997
  isl_stat (*add_ineq)(void *user, isl_int *), void *user)
1998
421
{
1999
421
  int i;
2000
421
  int r;
2001
421
  struct isl_mat *samples;
2002
421
  int nonneg;
2003
421
2004
421
  r = isl_tab_insert_div(tab, pos, div, add_ineq, user);
2005
421
  if (r < 0)
2006
0
    return isl_bool_error;
2007
421
  nonneg = tab->var[r].is_nonneg;
2008
421
  tab->var[r].frozen = 1;
2009
421
2010
421
  samples = isl_mat_extend(tab->samples,
2011
421
      tab->n_sample, 1 + tab->n_var);
2012
421
  tab->samples = samples;
2013
421
  if (!samples)
2014
0
    return isl_bool_error;
2015
1.64k
  
for (i = tab->n_outside; 421
i < samples->n_row1.64k
;
++i1.22k
)
{1.22k
2016
1.22k
    isl_seq_inner_product(div->el + 1, samples->row[i],
2017
1.22k
      div->size - 1, &samples->row[i][samples->n_col - 1]);
2018
1.22k
    isl_int_fdiv_q(samples->row[i][samples->n_col - 1],
2019
1.22k
             samples->row[i][samples->n_col - 1], div->el[0]);
2020
1.22k
  }
2021
421
  tab->samples = isl_mat_move_cols(tab->samples, 1 + pos,
2022
421
          1 + tab->n_var - 1, 1);
2023
421
  if (!tab->samples)
2024
0
    return isl_bool_error;
2025
421
2026
421
  return nonneg;
2027
421
}
2028
2029
/* Add a div specified by "div" to both the main tableau and
2030
 * the context tableau.  In case of the main tableau, we only
2031
 * need to add an extra div.  In the context tableau, we also
2032
 * need to express the meaning of the div.
2033
 * Return the index of the div or -1 if anything went wrong.
2034
 *
2035
 * The new integer division is added before any unknown integer
2036
 * divisions in the context to ensure that it does not get
2037
 * equated to some linear combination involving unknown integer
2038
 * divisions.
2039
 */
2040
static int add_div(struct isl_tab *tab, struct isl_context *context,
2041
  __isl_keep isl_vec *div)
2042
421
{
2043
421
  int r;
2044
421
  int pos;
2045
421
  isl_bool nonneg;
2046
421
  struct isl_tab *context_tab = context->op->peek_tab(context);
2047
421
2048
421
  if (
!tab || 421
!context_tab421
)
2049
0
    goto error;
2050
421
2051
421
  pos = context_tab->n_var - context->n_unknown;
2052
421
  if ((nonneg = context->op->insert_div(context, pos, div)) < 0)
2053
0
    goto error;
2054
421
2055
421
  
if (421
!context->op->is_ok(context)421
)
2056
0
    goto error;
2057
421
2058
421
  pos = tab->n_var - context->n_unknown;
2059
421
  if (isl_tab_extend_vars(tab, 1) < 0)
2060
0
    goto error;
2061
421
  r = isl_tab_insert_var(tab, pos);
2062
421
  if (r < 0)
2063
0
    goto error;
2064
421
  
if (421
nonneg421
)
2065
55
    tab->var[r].is_nonneg = 1;
2066
421
  tab->var[r].frozen = 1;
2067
421
  tab->n_div++;
2068
421
2069
421
  return tab->n_div - 1 - context->n_unknown;
2070
0
error:
2071
0
  context->op->invalidate(context);
2072
0
  return -1;
2073
421
}
2074
2075
static int find_div(struct isl_tab *tab, isl_int *div, isl_int denom)
2076
448
{
2077
448
  int i;
2078
448
  unsigned total = isl_basic_map_total_dim(tab->bmap);
2079
448
2080
1.64k
  for (i = 0; 
i < tab->bmap->n_div1.64k
;
++i1.19k
)
{1.22k
2081
1.22k
    if (isl_int_ne(tab->bmap->div[i][0], denom))
2082
863
      continue;
2083
363
    
if (363
!isl_seq_eq(tab->bmap->div[i] + 1, div, 1 + total)363
)
2084
336
      continue;
2085
27
    return i;
2086
363
  }
2087
421
  return -1;
2088
448
}
2089
2090
/* Return the index of a div that corresponds to "div".
2091
 * We first check if we already have such a div and if not, we create one.
2092
 */
2093
static int get_div(struct isl_tab *tab, struct isl_context *context,
2094
  struct isl_vec *div)
2095
448
{
2096
448
  int d;
2097
448
  struct isl_tab *context_tab = context->op->peek_tab(context);
2098
448
2099
448
  if (!context_tab)
2100
0
    return -1;
2101
448
2102
448
  d = find_div(context_tab, div->el + 1, div->el[0]);
2103
448
  if (d != -1)
2104
27
    return d;
2105
448
2106
421
  return add_div(tab, context, div);
2107
448
}
2108
2109
/* Add a parametric cut to cut away the non-integral sample value
2110
 * of the give row.
2111
 * Let a_i be the coefficients of the constant term and the parameters
2112
 * and let b_i be the coefficients of the variables or constraints
2113
 * in basis of the tableau.
2114
 * Let q be the div q = floor(\sum_i {-a_i} y_i).
2115
 *
2116
 * The cut is expressed as
2117
 *
2118
 *  c = \sum_i -{-a_i} y_i + \sum_i {b_i} x_i + q >= 0
2119
 *
2120
 * If q did not already exist in the context tableau, then it is added first.
2121
 * If q is in a column of the main tableau then the "+ q" can be accomplished
2122
 * by setting the corresponding entry to the denominator of the constraint.
2123
 * If q happens to be in a row of the main tableau, then the corresponding
2124
 * row needs to be added instead (taking care of the denominators).
2125
 * Note that this is very unlikely, but perhaps not entirely impossible.
2126
 *
2127
 * The current value of the cut is known to be negative (or at least
2128
 * non-positive), so row_sign is set accordingly.
2129
 *
2130
 * Return the row of the cut or -1.
2131
 */
2132
static int add_parametric_cut(struct isl_tab *tab, int row,
2133
  struct isl_context *context)
2134
372
{
2135
372
  struct isl_vec *div;
2136
372
  int d;
2137
372
  int i;
2138
372
  int r;
2139
372
  isl_int *r_row;
2140
372
  int col;
2141
372
  int n;
2142
372
  unsigned off = 2 + tab->M;
2143
372
2144
372
  if (!context)
2145
0
    return -1;
2146
372
2147
372
  div = get_row_parameter_div(tab, row);
2148
372
  if (!div)
2149
0
    return -1;
2150
372
2151
372
  n = tab->n_div - context->n_unknown;
2152
372
  d = context->op->get_div(context, tab, div);
2153
372
  isl_vec_free(div);
2154
372
  if (d < 0)
2155
0
    return -1;
2156
372
2157
372
  
if (372
isl_tab_extend_cons(tab, 1) < 0372
)
2158
0
    return -1;
2159
372
  r = isl_tab_allocate_con(tab);
2160
372
  if (r < 0)
2161
0
    return -1;
2162
372
2163
372
  r_row = tab->mat->row[tab->con[r].index];
2164
372
  isl_int_set(r_row[0], tab->mat->row[row][0]);
2165
372
  isl_int_neg(r_row[1], tab->mat->row[row][1]);
2166
372
  isl_int_fdiv_r(r_row[1], r_row[1], tab->mat->row[row][0]);
2167
372
  isl_int_neg(r_row[1], r_row[1]);
2168
372
  if (tab->M)
2169
372
    isl_int_set_si(r_row[2], 0);
2170
1.17k
  for (i = 0; 
i < tab->n_param1.17k
;
++i800
)
{800
2171
800
    if (tab->var[i].is_row)
2172
22
      continue;
2173
778
    col = tab->var[i].index;
2174
778
    isl_int_neg(r_row[off + col], tab->mat->row[row][off + col]);
2175
778
    isl_int_fdiv_r(r_row[off + col], r_row[off + col],
2176
778
        tab->mat->row[row][0]);
2177
778
    isl_int_neg(r_row[off + col], r_row[off + col]);
2178
778
  }
2179
1.88k
  for (i = 0; 
i < tab->n_div1.88k
;
++i1.51k
)
{1.51k
2180
1.51k
    if (tab->var[tab->n_var - tab->n_div + i].is_row)
2181
271
      continue;
2182
1.24k
    col = tab->var[tab->n_var - tab->n_div + i].index;
2183
1.24k
    isl_int_neg(r_row[off + col], tab->mat->row[row][off + col]);
2184
1.24k
    isl_int_fdiv_r(r_row[off + col], r_row[off + col],
2185
1.24k
        tab->mat->row[row][0]);
2186
1.24k
    isl_int_neg(r_row[off + col], r_row[off + col]);
2187
1.24k
  }
2188
3.81k
  for (i = 0; 
i < tab->n_col3.81k
;
++i3.44k
)
{3.44k
2189
3.44k
    if (tab->col_var[i] >= 0 &&
2190
2.01k
        (tab->col_var[i] < tab->n_param ||
2191
1.24k
         tab->col_var[i] >= tab->n_var - tab->n_div))
2192
2.01k
      continue;
2193
1.42k
    
isl_int_fdiv_r1.42k
(r_row[off + i],1.42k
2194
1.42k
      tab->mat->row[row][off + i], tab->mat->row[row][0]);
2195
1.42k
  }
2196
372
  if (
tab->var[tab->n_var - tab->n_div + d].is_row372
)
{1
2197
1
    isl_int gcd;
2198
1
    int d_row = tab->var[tab->n_var - tab->n_div + d].index;
2199
1
    isl_int_init(gcd);
2200
1
    isl_int_gcd(gcd, tab->mat->row[d_row][0], r_row[0]);
2201
1
    isl_int_divexact(r_row[0], r_row[0], gcd);
2202
1
    isl_int_divexact(gcd, tab->mat->row[d_row][0], gcd);
2203
1
    isl_seq_combine(r_row + 1, gcd, r_row + 1,
2204
1
        r_row[0], tab->mat->row[d_row] + 1,
2205
1
        off - 1 + tab->n_col);
2206
1
    isl_int_mul(r_row[0], r_row[0], tab->mat->row[d_row][0]);
2207
1
    isl_int_clear(gcd);
2208
371
  } else {
2209
371
    col = tab->var[tab->n_var - tab->n_div + d].index;
2210
371
    isl_int_set(r_row[off + col], tab->mat->row[row][0]);
2211
371
  }
2212
372
2213
372
  tab->con[r].is_nonneg = 1;
2214
372
  if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r]) < 0)
2215
0
    return -1;
2216
372
  
if (372
tab->row_sign372
)
2217
372
    tab->row_sign[tab->con[r].index] = isl_tab_row_neg;
2218
372
2219
372
  row = tab->con[r].index;
2220
372
2221
372
  if (
d >= n && 372
context->op->detect_equalities(context, tab) < 0370
)
2222
0
    return -1;
2223
372
2224
372
  return row;
2225
372
}
2226
2227
/* Construct a tableau for bmap that can be used for computing
2228
 * the lexicographic minimum (or maximum) of bmap.
2229
 * If not NULL, then dom is the domain where the minimum
2230
 * should be computed.  In this case, we set up a parametric
2231
 * tableau with row signs (initialized to "unknown").
2232
 * If M is set, then the tableau will use a big parameter.
2233
 * If max is set, then a maximum should be computed instead of a minimum.
2234
 * This means that for each variable x, the tableau will contain the variable
2235
 * x' = M - x, rather than x' = M + x.  This in turn means that the coefficient
2236
 * of the variables in all constraints are negated prior to adding them
2237
 * to the tableau.
2238
 */
2239
static __isl_give struct isl_tab *tab_for_lexmin(__isl_keep isl_basic_map *bmap,
2240
  __isl_keep isl_basic_set *dom, unsigned M, int max)
2241
4.21k
{
2242
4.21k
  int i;
2243
4.21k
  struct isl_tab *tab;
2244
4.21k
  unsigned n_var;
2245
4.21k
  unsigned o_var;
2246
4.21k
2247
4.21k
  tab = isl_tab_alloc(bmap->ctx, 2 * bmap->n_eq + bmap->n_ineq + 1,
2248
4.21k
          isl_basic_map_total_dim(bmap), M);
2249
4.21k
  if (!tab)
2250
0
    return NULL;
2251
4.21k
2252
4.21k
  
tab->rational = 4.21k
ISL_F_ISSET4.21k
(bmap, ISL_BASIC_MAP_RATIONAL);
2253
4.21k
  if (
dom4.21k
)
{3.80k
2254
3.80k
    tab->n_param = isl_basic_set_total_dim(dom) - dom->n_div;
2255
3.80k
    tab->n_div = dom->n_div;
2256
3.80k
    tab->row_sign = isl_calloc_array(bmap->ctx,
2257
3.80k
          enum isl_tab_row_sign, tab->mat->n_row);
2258
3.80k
    if (
tab->mat->n_row && 3.80k
!tab->row_sign3.80k
)
2259
0
      goto error;
2260
3.80k
  }
2261
4.21k
  
if (4.21k
ISL_F_ISSET4.21k
(bmap, ISL_BASIC_MAP_EMPTY))
{0
2262
0
    if (isl_tab_mark_empty(tab) < 0)
2263
0
      goto error;
2264
0
    return tab;
2265
0
  }
2266
4.21k
2267
18.8k
  
for (i = tab->n_param; 4.21k
i < tab->n_var - tab->n_div18.8k
;
++i14.6k
)
{14.6k
2268
14.6k
    tab->var[i].is_nonneg = 1;
2269
14.6k
    tab->var[i].frozen = 1;
2270
14.6k
  }
2271
4.21k
  o_var = 1 + tab->n_param;
2272
4.21k
  n_var = tab->n_var - tab->n_param - tab->n_div;
2273
19.3k
  for (i = 0; 
i < bmap->n_eq19.3k
;
++i15.1k
)
{15.1k
2274
15.1k
    if (max)
2275
12.3k
      isl_seq_neg(bmap->eq[i] + o_var,
2276
12.3k
            bmap->eq[i] + o_var, n_var);
2277
15.1k
    tab = add_lexmin_valid_eq(tab, bmap->eq[i]);
2278
15.1k
    if (max)
2279
12.3k
      isl_seq_neg(bmap->eq[i] + o_var,
2280
12.3k
            bmap->eq[i] + o_var, n_var);
2281
15.1k
    if (
!tab || 15.1k
tab->empty15.1k
)
2282
0
      return tab;
2283
15.1k
  }
2284
4.21k
  
if (4.21k
bmap->n_eq && 4.21k
restore_lexmin(tab) < 03.52k
)
2285
0
    goto error;
2286
17.9k
  
for (i = 0; 4.21k
i < bmap->n_ineq17.9k
;
++i13.7k
)
{13.7k
2287
13.7k
    if (max)
2288
8.36k
      isl_seq_neg(bmap->ineq[i] + o_var,
2289
8.36k
            bmap->ineq[i] + o_var, n_var);
2290
13.7k
    tab = add_lexmin_ineq(tab, bmap->ineq[i]);
2291
13.7k
    if (max)
2292
8.36k
      isl_seq_neg(bmap->ineq[i] + o_var,
2293
8.36k
            bmap->ineq[i] + o_var, n_var);
2294
13.7k
    if (
!tab || 13.7k
tab->empty13.7k
)
2295
0
      return tab;
2296
13.7k
  }
2297
4.21k
  return tab;
2298
0
error:
2299
0
  isl_tab_free(tab);
2300
0
  return NULL;
2301
4.21k
}
2302
2303
/* Given a main tableau where more than one row requires a split,
2304
 * determine and return the "best" row to split on.
2305
 *
2306
 * Given two rows in the main tableau, if the inequality corresponding
2307
 * to the first row is redundant with respect to that of the second row
2308
 * in the current tableau, then it is better to split on the second row,
2309
 * since in the positive part, both rows will be positive.
2310
 * (In the negative part a pivot will have to be performed and just about
2311
 * anything can happen to the sign of the other row.)
2312
 *
2313
 * As a simple heuristic, we therefore select the row that makes the most
2314
 * of the other rows redundant.
2315
 *
2316
 * Perhaps it would also be useful to look at the number of constraints
2317
 * that conflict with any given constraint.
2318
 *
2319
 * best is the best row so far (-1 when we have not found any row yet).
2320
 * best_r is the number of other rows made redundant by row best.
2321
 * When best is still -1, bset_r is meaningless, but it is initialized
2322
 * to some arbitrary value (0) anyway.  Without this redundant initialization
2323
 * valgrind may warn about uninitialized memory accesses when isl
2324
 * is compiled with some versions of gcc.
2325
 */
2326
static int best_split(struct isl_tab *tab, struct isl_tab *context_tab)
2327
303
{
2328
303
  struct isl_tab_undo *snap;
2329
303
  int split;
2330
303
  int row;
2331
303
  int best = -1;
2332
303
  int best_r = 0;
2333
303
2334
303
  if (isl_tab_extend_cons(context_tab, 2) < 0)
2335
0
    return -1;
2336
303
2337
303
  snap = isl_tab_snap(context_tab);
2338
303
2339
2.59k
  for (split = tab->n_redundant; 
split < tab->n_row2.59k
;
++split2.28k
)
{2.28k
2340
2.28k
    struct isl_tab_undo *snap2;
2341
2.28k
    struct isl_vec *ineq = NULL;
2342
2.28k
    int r = 0;
2343
2.28k
    int ok;
2344
2.28k
2345
2.28k
    if (!isl_tab_var_from_row(tab, split)->is_nonneg)
2346
34
      continue;
2347
2.25k
    
if (2.25k
tab->row_sign[split] != isl_tab_row_any2.25k
)
2348
1.54k
      continue;
2349
2.25k
2350
709
    ineq = get_row_parameter_ineq(tab, split);
2351
709
    if (!ineq)
2352
0
      return -1;
2353
709
    ok = isl_tab_add_ineq(context_tab, ineq->el) >= 0;
2354
709
    isl_vec_free(ineq);
2355
709
    if (!ok)
2356
0
      return -1;
2357
709
2358
709
    snap2 = isl_tab_snap(context_tab);
2359
709
2360
6.31k
    for (row = tab->n_redundant; 
row < tab->n_row6.31k
;
++row5.60k
)
{5.60k
2361
5.60k
      struct isl_tab_var *var;
2362
5.60k
2363
5.60k
      if (row == split)
2364
709
        continue;
2365
4.89k
      
if (4.89k
!isl_tab_var_from_row(tab, row)->is_nonneg4.89k
)
2366
74
        continue;
2367
4.81k
      
if (4.81k
tab->row_sign[row] != isl_tab_row_any4.81k
)
2368
3.64k
        continue;
2369
4.81k
2370
1.17k
      ineq = get_row_parameter_ineq(tab, row);
2371
1.17k
      if (!ineq)
2372
0
        return -1;
2373
1.17k
      ok = isl_tab_add_ineq(context_tab, ineq->el) >= 0;
2374
1.17k
      isl_vec_free(ineq);
2375
1.17k
      if (!ok)
2376
0
        return -1;
2377
1.17k
      var = &context_tab->con[context_tab->n_con - 1];
2378
1.17k
      if (!context_tab->empty &&
2379
1.17k
          !isl_tab_min_at_most_neg_one(context_tab, var))
2380
51
        r++;
2381
1.17k
      if (isl_tab_rollback(context_tab, snap2) < 0)
2382
0
        return -1;
2383
1.17k
    }
2384
709
    
if (709
best == -1 || 709
r > best_r406
)
{321
2385
321
      best = split;
2386
321
      best_r = r;
2387
321
    }
2388
709
    if (isl_tab_rollback(context_tab, snap) < 0)
2389
0
      return -1;
2390
709
  }
2391
303
2392
303
  return best;
2393
303
}
2394
2395
static struct isl_basic_set *context_lex_peek_basic_set(
2396
  struct isl_context *context)
2397
0
{
2398
0
  struct isl_context_lex *clex = (struct isl_context_lex *)context;
2399
0
  if (!clex->tab)
2400
0
    return NULL;
2401
0
  return isl_tab_peek_bset(clex->tab);
2402
0
}
2403
2404
static struct isl_tab *context_lex_peek_tab(struct isl_context *context)
2405
0
{
2406
0
  struct isl_context_lex *clex = (struct isl_context_lex *)context;
2407
0
  return clex->tab;
2408
0
}
2409
2410
static void context_lex_add_eq(struct isl_context *context, isl_int *eq,
2411
    int check, int update)
2412
0
{
2413
0
  struct isl_context_lex *clex = (struct isl_context_lex *)context;
2414
0
  if (isl_tab_extend_cons(clex->tab, 2) < 0)
2415
0
    goto error;
2416
0
  
if (0
add_lexmin_eq(clex->tab, eq) < 00
)
2417
0
    goto error;
2418
0
  
if (0
check0
)
{0
2419
0
    int v = tab_has_valid_sample(clex->tab, eq, 1);
2420
0
    if (v < 0)
2421
0
      goto error;
2422
0
    
if (0
!v0
)
2423
0
      clex->tab = check_integer_feasible(clex->tab);
2424
0
  }
2425
0
  
if (0
update0
)
2426
0
    clex->tab = check_samples(clex->tab, eq, 1);
2427
0
  return;
2428
0
error:
2429
0
  isl_tab_free(clex->tab);
2430
0
  clex->tab = NULL;
2431
0
}
2432
2433
static void context_lex_add_ineq(struct isl_context *context, isl_int *ineq,
2434
    int check, int update)
2435
0
{
2436
0
  struct isl_context_lex *clex = (struct isl_context_lex *)context;
2437
0
  if (isl_tab_extend_cons(clex->tab, 1) < 0)
2438
0
    goto error;
2439
0
  clex->tab = add_lexmin_ineq(clex->tab, ineq);
2440
0
  if (
check0
)
{0
2441
0
    int v = tab_has_valid_sample(clex->tab, ineq, 0);
2442
0
    if (v < 0)
2443
0
      goto error;
2444
0
    
if (0
!v0
)
2445
0
      clex->tab = check_integer_feasible(clex->tab);
2446
0
  }
2447
0
  
if (0
update0
)
2448
0
    clex->tab = check_samples(clex->tab, ineq, 0);
2449
0
  return;
2450
0
error:
2451
0
  isl_tab_free(clex->tab);
2452
0
  clex->tab = NULL;
2453
0
}
2454
2455
static isl_stat context_lex_add_ineq_wrap(void *user, isl_int *ineq)
2456
0
{
2457
0
  struct isl_context *context = (struct isl_context *)user;
2458
0
  context_lex_add_ineq(context, ineq, 0, 0);
2459
0
  return context->op->is_ok(context) ? 
isl_stat_ok0
:
isl_stat_error0
;
2460
0
}
2461
2462
/* Check which signs can be obtained by "ineq" on all the currently
2463
 * active sample values.  See row_sign for more information.
2464
 */
2465
static enum isl_tab_row_sign tab_ineq_sign(struct isl_tab *tab, isl_int *ineq,
2466
  int strict)
2467
7.82k
{
2468
7.82k
  int i;
2469
7.82k
  int sgn;
2470
7.82k
  isl_int tmp;
2471
7.82k
  enum isl_tab_row_sign res = isl_tab_row_unknown;
2472
7.82k
2473
7.82k
  isl_assert(tab->mat->ctx, tab->samples, return isl_tab_row_unknown);
2474
7.82k
  
isl_assert7.82k
(tab->mat->ctx, tab->samples->n_col == 1 + tab->n_var,7.82k
2475
7.82k
      return isl_tab_row_unknown);
2476
7.82k
2477
7.82k
  
isl_int_init7.82k
(tmp);7.82k
2478
19.3k
  for (i = tab->n_outside; 
i < tab->n_sample19.3k
;
++i11.4k
)
{11.9k
2479
11.9k
    isl_seq_inner_product(tab->samples->row[i], ineq,
2480
11.9k
          1 + tab->n_var, &tmp);
2481
11.9k
    sgn = isl_int_sgn(tmp);
2482
11.9k
    if (
sgn > 0 || 11.9k
(sgn == 0 && 5.33k
strict2.91k
))
{9.32k
2483
9.32k
      if (res == isl_tab_row_unknown)
2484
5.94k
        res = isl_tab_row_pos;
2485
9.32k
      if (res == isl_tab_row_neg)
2486
214
        res = isl_tab_row_any;
2487
9.32k
    }
2488
11.9k
    if (
sgn < 011.9k
)
{2.42k
2489
2.42k
      if (res == isl_tab_row_unknown)
2490
1.84k
        res = isl_tab_row_neg;
2491
2.42k
      if (res == isl_tab_row_pos)
2492
240
        res = isl_tab_row_any;
2493
2.42k
    }
2494
11.9k
    if (res == isl_tab_row_any)
2495
454
      break;
2496
11.9k
  }
2497
7.82k
  isl_int_clear(tmp);
2498
7.82k
2499
7.82k
  return res;
2500
7.82k
}
2501
2502
static enum isl_tab_row_sign context_lex_ineq_sign(struct isl_context *context,
2503
      isl_int *ineq, int strict)
2504
0
{
2505
0
  struct isl_context_lex *clex = (struct isl_context_lex *)context;
2506
0
  return tab_ineq_sign(clex->tab, ineq, strict);
2507
0
}
2508
2509
/* Check whether "ineq" can be added to the tableau without rendering
2510
 * it infeasible.
2511
 */
2512
static int context_lex_test_ineq(struct isl_context *context, isl_int *ineq)
2513
0
{
2514
0
  struct isl_context_lex *clex = (struct isl_context_lex *)context;
2515
0
  struct isl_tab_undo *snap;
2516
0
  int feasible;
2517
0
2518
0
  if (!clex->tab)
2519
0
    return -1;
2520
0
2521
0
  
if (0
isl_tab_extend_cons(clex->tab, 1) < 00
)
2522
0
    return -1;
2523
0
2524
0
  snap = isl_tab_snap(clex->tab);
2525
0
  if (isl_tab_push_basis(clex->tab) < 0)
2526
0
    return -1;
2527
0
  clex->tab = add_lexmin_ineq(clex->tab, ineq);
2528
0
  clex->tab = check_integer_feasible(clex->tab);
2529
0
  if (!clex->tab)
2530
0
    return -1;
2531
0
  feasible = !clex->tab->empty;
2532
0
  if (isl_tab_rollback(clex->tab, snap) < 0)
2533
0
    return -1;
2534
0
2535
0
  return feasible;
2536
0
}
2537
2538
static int context_lex_get_div(struct isl_context *context, struct isl_tab *tab,
2539
    struct isl_vec *div)
2540
0
{
2541
0
  return get_div(tab, context, div);
2542
0
}
2543
2544
/* Insert a div specified by "div" to the context tableau at position "pos" and
2545
 * return isl_bool_true if the div is obviously non-negative.
2546
 * context_tab_add_div will always return isl_bool_true, because all variables
2547
 * in a isl_context_lex tableau are non-negative.
2548
 * However, if we are using a big parameter in the context, then this only
2549
 * reflects the non-negativity of the variable used to _encode_ the
2550
 * div, i.e., div' = M + div, so we can't draw any conclusions.
2551
 */
2552
static isl_bool context_lex_insert_div(struct isl_context *context, int pos,
2553
  __isl_keep isl_vec *div)
2554
0
{
2555
0
  struct isl_context_lex *clex = (struct isl_context_lex *)context;
2556
0
  isl_bool nonneg;
2557
0
  nonneg = context_tab_insert_div(clex->tab, pos, div,
2558
0
          context_lex_add_ineq_wrap, context);
2559
0
  if (nonneg < 0)
2560
0
    return isl_bool_error;
2561
0
  
if (0
clex->tab->M0
)
2562
0
    return isl_bool_false;
2563
0
  return nonneg;
2564
0
}
2565
2566
static int context_lex_detect_equalities(struct isl_context *context,
2567
    struct isl_tab *tab)
2568
0
{
2569
0
  return 0;
2570
0
}
2571
2572
static int context_lex_best_split(struct isl_context *context,
2573
    struct isl_tab *tab)
2574
0
{
2575
0
  struct isl_context_lex *clex = (struct isl_context_lex *)context;
2576
0
  struct isl_tab_undo *snap;
2577
0
  int r;
2578
0
2579
0
  snap = isl_tab_snap(clex->tab);
2580
0
  if (isl_tab_push_basis(clex->tab) < 0)
2581
0
    return -1;
2582
0
  r = best_split(tab, clex->tab);
2583
0
2584
0
  if (
r >= 0 && 0
isl_tab_rollback(clex->tab, snap) < 00
)
2585
0
    return -1;
2586
0
2587
0
  return r;
2588
0
}
2589
2590
static int context_lex_is_empty(struct isl_context *context)
2591
0
{
2592
0
  struct isl_context_lex *clex = (struct isl_context_lex *)context;
2593
0
  if (!clex->tab)
2594
0
    return -1;
2595
0
  return clex->tab->empty;
2596
0
}
2597
2598
static void *context_lex_save(struct isl_context *context)
2599
0
{
2600
0
  struct isl_context_lex *clex = (struct isl_context_lex *)context;
2601
0
  struct isl_tab_undo *snap;
2602
0
2603
0
  snap = isl_tab_snap(clex->tab);
2604
0
  if (isl_tab_push_basis(clex->tab) < 0)
2605
0
    return NULL;
2606
0
  
if (0
isl_tab_save_samples(clex->tab) < 00
)
2607
0
    return NULL;
2608
0
2609
0
  return snap;
2610
0
}
2611
2612
static void context_lex_restore(struct isl_context *context, void *save)
2613
0
{
2614
0
  struct isl_context_lex *clex = (struct isl_context_lex *)context;
2615
0
  if (
isl_tab_rollback(clex->tab, (struct isl_tab_undo *)save) < 00
)
{0
2616
0
    isl_tab_free(clex->tab);
2617
0
    clex->tab = NULL;
2618
0
  }
2619
0
}
2620
2621
static void context_lex_discard(void *save)
2622
0
{
2623
0
}
2624
2625
static int context_lex_is_ok(struct isl_context *context)
2626
0
{
2627
0
  struct isl_context_lex *clex = (struct isl_context_lex *)context;
2628
0
  return !!clex->tab;
2629
0
}
2630
2631
/* For each variable in the context tableau, check if the variable can
2632
 * only attain non-negative values.  If so, mark the parameter as non-negative
2633
 * in the main tableau.  This allows for a more direct identification of some
2634
 * cases of violated constraints.
2635
 */
2636
static struct isl_tab *tab_detect_nonnegative_parameters(struct isl_tab *tab,
2637
  struct isl_tab *context_tab)
2638
3.80k
{
2639
3.80k
  int i;
2640
3.80k
  struct isl_tab_undo *snap;
2641
3.80k
  struct isl_vec *ineq = NULL;
2642
3.80k
  struct isl_tab_var *var;
2643
3.80k
  int n;
2644
3.80k
2645
3.80k
  if (context_tab->n_var == 0)
2646
609
    return tab;
2647
3.80k
2648
3.19k
  ineq = isl_vec_alloc(tab->mat->ctx, 1 + context_tab->n_var);
2649
3.19k
  if (!ineq)
2650
0
    goto error;
2651
3.19k
2652
3.19k
  
if (3.19k
isl_tab_extend_cons(context_tab, 1) < 03.19k
)
2653
0
    goto error;
2654
3.19k
2655
3.19k
  snap = isl_tab_snap(context_tab);
2656
3.19k
2657
3.19k
  n = 0;
2658
3.19k
  isl_seq_clr(ineq->el, ineq->size);
2659
16.5k
  for (i = 0; 
i < context_tab->n_var16.5k
;
++i13.3k
)
{13.3k
2660
13.3k
    isl_int_set_si(ineq->el[1 + i], 1);
2661
13.3k
    if (isl_tab_add_ineq(context_tab, ineq->el) < 0)
2662
0
      goto error;
2663
13.3k
    var = &context_tab->con[context_tab->n_con - 1];
2664
13.3k
    if (!context_tab->empty &&
2665
13.1k
        
!isl_tab_min_at_most_neg_one(context_tab, var)13.1k
)
{10.0k
2666
10.0k
      int j = i;
2667
10.0k
      if (i >= tab->n_param)
2668
41
        j = i - tab->n_param + tab->n_var - tab->n_div;
2669
10.0k
      tab->var[j].is_nonneg = 1;
2670
10.0k
      n++;
2671
10.0k
    }
2672
13.3k
    isl_int_set_si(ineq->el[1 + i], 0);
2673
13.3k
    if (isl_tab_rollback(context_tab, snap) < 0)
2674
0
      goto error;
2675
13.3k
  }
2676
3.19k
2677
3.19k
  
if (3.19k
context_tab->M && 3.19k
n == context_tab->n_var0
)
{0
2678
0
    context_tab->mat = isl_mat_drop_cols(context_tab->mat, 2, 1);
2679
0
    context_tab->M = 0;
2680
0
  }
2681
3.19k
2682
3.19k
  isl_vec_free(ineq);
2683
3.19k
  return tab;
2684
0
error:
2685
0
  isl_vec_free(ineq);
2686
0
  isl_tab_free(tab);
2687
0
  return NULL;
2688
3.19k
}
2689
2690
static struct isl_tab *context_lex_detect_nonnegative_parameters(
2691
  struct isl_context *context, struct isl_tab *tab)
2692
0
{
2693
0
  struct isl_context_lex *clex = (struct isl_context_lex *)context;
2694
0
  struct isl_tab_undo *snap;
2695
0
2696
0
  if (!tab)
2697
0
    return NULL;
2698
0
2699
0
  snap = isl_tab_snap(clex->tab);
2700
0
  if (isl_tab_push_basis(clex->tab) < 0)
2701
0
    goto error;
2702
0
2703
0
  tab = tab_detect_nonnegative_parameters(tab, clex->tab);
2704
0
2705
0
  if (isl_tab_rollback(clex->tab, snap) < 0)
2706
0
    goto error;
2707
0
2708
0
  return tab;
2709
0
error:
2710
0
  isl_tab_free(tab);
2711
0
  return NULL;
2712
0
}
2713
2714
static void context_lex_invalidate(struct isl_context *context)
2715
0
{
2716
0
  struct isl_context_lex *clex = (struct isl_context_lex *)context;
2717
0
  isl_tab_free(clex->tab);
2718
0
  clex->tab = NULL;
2719
0
}
2720
2721
static __isl_null struct isl_context *context_lex_free(
2722
  struct isl_context *context)
2723
0
{
2724
0
  struct isl_context_lex *clex = (struct isl_context_lex *)context;
2725
0
  isl_tab_free(clex->tab);
2726
0
  free(clex);
2727
0
2728
0
  return NULL;
2729
0
}
2730
2731
struct isl_context_op isl_context_lex_op = {
2732
  context_lex_detect_nonnegative_parameters,
2733
  context_lex_peek_basic_set,
2734
  context_lex_peek_tab,
2735
  context_lex_add_eq,
2736
  context_lex_add_ineq,
2737
  context_lex_ineq_sign,
2738
  context_lex_test_ineq,
2739
  context_lex_get_div,
2740
  context_lex_insert_div,
2741
  context_lex_detect_equalities,
2742
  context_lex_best_split,
2743
  context_lex_is_empty,
2744
  context_lex_is_ok,
2745
  context_lex_save,
2746
  context_lex_restore,
2747
  context_lex_discard,
2748
  context_lex_invalidate,
2749
  context_lex_free,
2750
};
2751
2752
static struct isl_tab *context_tab_for_lexmin(__isl_take isl_basic_set *bset)
2753
0
{
2754
0
  struct isl_tab *tab;
2755
0
2756
0
  if (!bset)
2757
0
    return NULL;
2758
0
  tab = tab_for_lexmin(bset_to_bmap(bset), NULL, 1, 0);
2759
0
  if (isl_tab_track_bset(tab, bset) < 0)
2760
0
    goto error;
2761
0
  tab = isl_tab_init_samples(tab);
2762
0
  return tab;
2763
0
error:
2764
0
  isl_tab_free(tab);
2765
0
  return NULL;
2766
0
}
2767
2768
static struct isl_context *isl_context_lex_alloc(struct isl_basic_set *dom)
2769
0
{
2770
0
  struct isl_context_lex *clex;
2771
0
2772
0
  if (!dom)
2773
0
    return NULL;
2774
0
2775
0
  
clex = 0
isl_alloc_type0
(dom->ctx, struct isl_context_lex);
2776
0
  if (!clex)
2777
0
    return NULL;
2778
0
2779
0
  clex->context.op = &isl_context_lex_op;
2780
0
2781
0
  clex->tab = context_tab_for_lexmin(isl_basic_set_copy(dom));
2782
0
  if (restore_lexmin(clex->tab) < 0)
2783
0
    goto error;
2784
0
  clex->tab = check_integer_feasible(clex->tab);
2785
0
  if (!clex->tab)
2786
0
    goto error;
2787
0
2788
0
  return &clex->context;
2789
0
error:
2790
0
  clex->context.op->free(&clex->context);
2791
0
  return NULL;
2792
0
}
2793
2794
/* Representation of the context when using generalized basis reduction.
2795
 *
2796
 * "shifted" contains the offsets of the unit hypercubes that lie inside the
2797
 * context.  Any rational point in "shifted" can therefore be rounded
2798
 * up to an integer point in the context.
2799
 * If the context is constrained by any equality, then "shifted" is not used
2800
 * as it would be empty.
2801
 */
2802
struct isl_context_gbr {
2803
  struct isl_context context;
2804
  struct isl_tab *tab;
2805
  struct isl_tab *shifted;
2806
  struct isl_tab *cone;
2807
};
2808
2809
static struct isl_tab *context_gbr_detect_nonnegative_parameters(
2810
  struct isl_context *context, struct isl_tab *tab)
2811
3.80k
{
2812
3.80k
  struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
2813
3.80k
  if (!tab)
2814
0
    return NULL;
2815
3.80k
  return tab_detect_nonnegative_parameters(tab, cgbr->tab);
2816
3.80k
}
2817
2818
static struct isl_basic_set *context_gbr_peek_basic_set(
2819
  struct isl_context *context)
2820
15.3k
{
2821
15.3k
  struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
2822
15.3k
  if (!cgbr->tab)
2823
0
    return NULL;
2824
15.3k
  return isl_tab_peek_bset(cgbr->tab);
2825
15.3k
}
2826
2827
static struct isl_tab *context_gbr_peek_tab(struct isl_context *context)
2828
21.9k
{
2829
21.9k
  struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
2830
21.9k
  return cgbr->tab;
2831
21.9k
}
2832
2833
/* Initialize the "shifted" tableau of the context, which
2834
 * contains the constraints of the original tableau shifted
2835
 * by the sum of all negative coefficients.  This ensures
2836
 * that any rational point in the shifted tableau can
2837
 * be rounded up to yield an integer point in the original tableau.
2838
 */
2839
static void gbr_init_shifted(struct isl_context_gbr *cgbr)
2840
49
{
2841
49
  int i, j;
2842
49
  struct isl_vec *cst;
2843
49
  struct isl_basic_set *bset = isl_tab_peek_bset(cgbr->tab);
2844
49
  unsigned dim = isl_basic_set_total_dim(bset);
2845
49
2846
49
  cst = isl_vec_alloc(cgbr->tab->mat->ctx, bset->n_ineq);
2847
49
  if (!cst)
2848
0
    return;
2849
49
2850
576
  
for (i = 0; 49
i < bset->n_ineq576
;
++i527
)
{527
2851
527
    isl_int_set(cst->el[i], bset->ineq[i][0]);
2852
5.42k
    for (j = 0; 
j < dim5.42k
;
++j4.89k
)
{4.89k
2853
4.89k
      if (
!4.89k
isl_int_is_neg4.89k
(bset->ineq[i][1 + j]))
2854
4.42k
        continue;
2855
476
      
isl_int_add476
(bset->ineq[i][0], bset->ineq[i][0],476
2856
476
            bset->ineq[i][1 + j]);
2857
476
    }
2858
527
  }
2859
49
2860
49
  cgbr->shifted = isl_tab_from_basic_set(bset, 0);
2861
49
2862
576
  for (i = 0; 
i < bset->n_ineq576
;
++i527
)
2863
527
    isl_int_set(bset->ineq[i][0], cst->el[i]);
2864
49
2865
49
  isl_vec_free(cst);
2866
49
}
2867
2868
/* Check if the shifted tableau is non-empty, and if so
2869
 * use the sample point to construct an integer point
2870
 * of the context tableau.
2871
 */
2872
static struct isl_vec *gbr_get_shifted_sample(struct isl_context_gbr *cgbr)
2873
54
{
2874
54
  struct isl_vec *sample;
2875
54
2876
54
  if (!cgbr->shifted)
2877
49
    gbr_init_shifted(cgbr);
2878
54
  if (!cgbr->shifted)
2879
0
    return NULL;
2880
54
  
if (54
cgbr->shifted->empty54
)
2881
34
    return isl_vec_alloc(cgbr->tab->mat->ctx, 0);
2882
54
2883
20
  sample = isl_tab_get_sample_value(cgbr->shifted);
2884
20
  sample = isl_vec_ceil(sample);
2885
20
2886
20
  return sample;
2887
54
}
2888
2889
static __isl_give isl_basic_set *drop_constant_terms(
2890
  __isl_take isl_basic_set *bset)
2891
607
{
2892
607
  int i;
2893
607
2894
607
  if (!bset)
2895
0
    return NULL;
2896
607
2897
1.66k
  
for (i = 0; 607
i < bset->n_eq1.66k
;
++i1.05k
)
2898
1.05k
    isl_int_set_si(bset->eq[i][0], 0);
2899
607
2900
4.36k
  for (i = 0; 
i < bset->n_ineq4.36k
;
++i3.75k
)
2901
3.75k
    isl_int_set_si(bset->ineq[i][0], 0);
2902
607
2903
607
  return bset;
2904
607
}
2905
2906
static int use_shifted(struct isl_context_gbr *cgbr)
2907
708
{
2908
708
  if (!cgbr->tab)
2909
0
    return 0;
2910
708
  
return cgbr->tab->bmap->n_eq == 0 && 708
cgbr->tab->bmap->n_div == 0154
;
2911
708
}
2912
2913
static struct isl_vec *gbr_get_sample(struct isl_context_gbr *cgbr)
2914
7.10k
{
2915
7.10k
  struct isl_basic_set *bset;
2916
7.10k
  struct isl_basic_set *cone;
2917
7.10k
2918
7.10k
  if (isl_tab_sample_is_integer(cgbr->tab))
2919
6.40k
    return isl_tab_get_sample_value(cgbr->tab);
2920
7.10k
2921
705
  
if (705
use_shifted(cgbr)705
)
{54
2922
54
    struct isl_vec *sample;
2923
54
2924
54
    sample = gbr_get_shifted_sample(cgbr);
2925
54
    if (
!sample || 54
sample->size > 054
)
2926
20
      return sample;
2927
54
2928
34
    isl_vec_free(sample);
2929
34
  }
2930
705
2931
685
  
if (685
!cgbr->cone685
)
{556
2932
556
    bset = isl_tab_peek_bset(cgbr->tab);
2933
556
    cgbr->cone = isl_tab_from_recession_cone(bset, 0);
2934
556
    if (!cgbr->cone)
2935
0
      return NULL;
2936
556
    
if (556
isl_tab_track_bset(cgbr->cone,556
2937
556
          isl_basic_set_copy(bset)) < 0)
2938
0
      return NULL;
2939
556
  }
2940
685
  
if (685
isl_tab_detect_implicit_equalities(cgbr->cone) < 0685
)
2941
0
    return NULL;
2942
685
2943
685
  
if (685
cgbr->cone->n_dead == cgbr->cone->n_col685
)
{78
2944
78
    struct isl_vec *sample;
2945
78
    struct isl_tab_undo *snap;
2946
78
2947
78
    if (
cgbr->tab->basis78
)
{30
2948
30
      if (
cgbr->tab->basis->n_col != 1 + cgbr->tab->n_var30
)
{12
2949
12
        isl_mat_free(cgbr->tab->basis);
2950
12
        cgbr->tab->basis = NULL;
2951
12
      }
2952
30
      cgbr->tab->n_zero = 0;
2953
30
      cgbr->tab->n_unbounded = 0;
2954
30
    }
2955
78
2956
78
    snap = isl_tab_snap(cgbr->tab);
2957
78
2958
78
    sample = isl_tab_sample(cgbr->tab);
2959
78
2960
78
    if (
!sample || 78
isl_tab_rollback(cgbr->tab, snap) < 078
)
{0
2961
0
      isl_vec_free(sample);
2962
0
      return NULL;
2963
0
    }
2964
78
2965
78
    return sample;
2966
78
  }
2967
685
2968
607
  cone = isl_basic_set_dup(isl_tab_peek_bset(cgbr->cone));
2969
607
  cone = drop_constant_terms(cone);
2970
607
  cone = isl_basic_set_update_from_tab(cone, cgbr->cone);
2971
607
  cone = isl_basic_set_underlying_set(cone);
2972
607
  cone = isl_basic_set_gauss(cone, NULL);
2973
607
2974
607
  bset = isl_basic_set_dup(isl_tab_peek_bset(cgbr->tab));
2975
607
  bset = isl_basic_set_update_from_tab(bset, cgbr->tab);
2976
607
  bset = isl_basic_set_underlying_set(bset);
2977
607
  bset = isl_basic_set_gauss(bset, NULL);
2978
607
2979
607
  return isl_basic_set_sample_with_cone(bset, cone);
2980
685
}
2981
2982
static void check_gbr_integer_feasible(struct isl_context_gbr *cgbr)
2983
25.4k
{
2984
25.4k
  struct isl_vec *sample;
2985
25.4k
2986
25.4k
  if (!cgbr->tab)
2987
0
    return;
2988
25.4k
2989
25.4k
  
if (25.4k
cgbr->tab->empty25.4k
)
2990
18.3k
    return;
2991
25.4k
2992
7.10k
  sample = gbr_get_sample(cgbr);
2993
7.10k
  if (!sample)
2994
0
    goto error;
2995
7.10k
2996
7.10k
  
if (7.10k
sample->size == 07.10k
)
{132
2997
132
    isl_vec_free(sample);
2998
132
    if (isl_tab_mark_empty(cgbr->tab) < 0)
2999
0
      goto error;
3000
132
    return;
3001
132
  }
3002
7.10k
3003
6.97k
  
if (6.97k
isl_tab_add_sample(cgbr->tab, sample) < 06.97k
)
3004
0
    goto error;
3005
6.97k
3006
6.97k
  return;
3007
0
error:
3008
0
  isl_tab_free(cgbr->tab);
3009
0
  cgbr->tab = NULL;
3010
0
}
3011
3012
static struct isl_tab *add_gbr_eq(struct isl_tab *tab, isl_int *eq)
3013
6.60k
{
3014
6.60k
  if (!tab)
3015
0
    return NULL;
3016
6.60k
3017
6.60k
  
if (6.60k
isl_tab_extend_cons(tab, 2) < 06.60k
)
3018
0
    goto error;
3019
6.60k
3020
6.60k
  
if (6.60k
isl_tab_add_eq(tab, eq) < 06.60k
)
3021
0
    goto error;
3022
6.60k
3023
6.60k
  return tab;
3024
0
error:
3025
0
  isl_tab_free(tab);
3026
0
  return NULL;
3027
6.60k
}
3028
3029
/* Add the equality described by "eq" to the context.
3030
 * If "check" is set, then we check if the context is empty after
3031
 * adding the equality.
3032
 * If "update" is set, then we check if the samples are still valid.
3033
 *
3034
 * We do not explicitly add shifted copies of the equality to
3035
 * cgbr->shifted since they would conflict with each other.
3036
 * Instead, we directly mark cgbr->shifted empty.
3037
 */
3038
static void context_gbr_add_eq(struct isl_context *context, isl_int *eq,
3039
    int check, int update)
3040
6.60k
{
3041
6.60k
  struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3042
6.60k
3043
6.60k
  cgbr->tab = add_gbr_eq(cgbr->tab, eq);
3044
6.60k
3045
6.60k
  if (
cgbr->shifted && 6.60k
!cgbr->shifted->empty2
&&
use_shifted(cgbr)1
)
{1
3046
1
    if (isl_tab_mark_empty(cgbr->shifted) < 0)
3047
0
      goto error;
3048
1
  }
3049
6.60k
3050
6.60k
  
if (6.60k
cgbr->cone && 6.60k
cgbr->cone->n_col != cgbr->cone->n_dead221
)
{206
3051
206
    if (isl_tab_extend_cons(cgbr->cone, 2) < 0)
3052
0
      goto error;
3053
206
    
if (206
isl_tab_add_eq(cgbr->cone, eq) < 0206
)
3054
0
      goto error;
3055
206
  }
3056
6.60k
3057
6.60k
  
if (6.60k
check6.60k
)
{6.60k
3058
6.60k
    int v = tab_has_valid_sample(cgbr->tab, eq, 1);
3059
6.60k
    if (v < 0)
3060
0
      goto error;
3061
6.60k
    
if (6.60k
!v6.60k
)
3062
83
      check_gbr_integer_feasible(cgbr);
3063
6.60k
  }
3064
6.60k
  
if (6.60k
update6.60k
)
3065
6.60k
    cgbr->tab = check_samples(cgbr->tab, eq, 1);
3066
6.60k
  return;
3067
0
error:
3068
0
  isl_tab_free(cgbr->tab);
3069
0
  cgbr->tab = NULL;
3070
0
}
3071
3072
static void add_gbr_ineq(struct isl_context_gbr *cgbr, isl_int *ineq)
3073
25.4k
{
3074
25.4k
  if (!cgbr->tab)
3075
0
    return;
3076
25.4k
3077
25.4k
  
if (25.4k
isl_tab_extend_cons(cgbr->tab, 1) < 025.4k
)
3078
0
    goto error;
3079
25.4k
3080
25.4k
  
if (25.4k
isl_tab_add_ineq(cgbr->tab, ineq) < 025.4k
)
3081
0
    goto error;
3082
25.4k
3083
25.4k
  
if (25.4k
cgbr->shifted && 25.4k
!cgbr->shifted->empty30
&&
use_shifted(cgbr)2
)
{2
3084
2
    int i;
3085
2
    unsigned dim;
3086
2
    dim = isl_basic_map_total_dim(cgbr->tab->bmap);
3087
2
3088
2
    if (isl_tab_extend_cons(cgbr->shifted, 1) < 0)
3089
0
      goto error;
3090
2
3091
10
    
for (i = 0; 2
i < dim10
;
++i8
)
{8
3092
8
      if (
!8
isl_int_is_neg8
(ineq[1 + i]))
3093
6
        continue;
3094
2
      
isl_int_add2
(ineq[0], ineq[0], ineq[1 + i]);2
3095
2
    }
3096
2
3097
2
    if (isl_tab_add_ineq(cgbr->shifted, ineq) < 0)
3098
0
      goto error;
3099
2
3100
10
    
for (i = 0; 2
i < dim10
;
++i8
)
{8
3101
8
      if (
!8
isl_int_is_neg8
(ineq[1 + i]))
3102
6
        continue;
3103
2
      
isl_int_sub2
(ineq[0], ineq[0], ineq[1 + i]);2
3104
2
    }
3105
2
  }
3106
25.4k
3107
25.4k
  
if (25.4k
cgbr->cone && 25.4k
cgbr->cone->n_col != cgbr->cone->n_dead3.19k
)
{2.17k
3108
2.17k
    if (isl_tab_extend_cons(cgbr->cone, 1) < 0)
3109
0
      goto error;
3110
2.17k
    
if (2.17k
isl_tab_add_ineq(cgbr->cone, ineq) < 02.17k
)
3111
0
      goto error;
3112
2.17k
  }
3113
25.4k
3114
25.4k
  return;
3115
0
error:
3116
0
  isl_tab_free(cgbr->tab);
3117
0
  cgbr->tab = NULL;
3118
0
}
3119
3120
static void context_gbr_add_ineq(struct isl_context *context, isl_int *ineq,
3121
    int check, int update)
3122
18.1k
{
3123
18.1k
  struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3124
18.1k
3125
18.1k
  add_gbr_ineq(cgbr, ineq);
3126
18.1k
  if (!cgbr->tab)
3127
0
    return;
3128
18.1k
3129
18.1k
  
if (18.1k
check18.1k
)
{13.3k
3130
13.3k
    int v = tab_has_valid_sample(cgbr->tab, ineq, 0);
3131
13.3k
    if (v < 0)
3132
0
      goto error;
3133
13.3k
    
if (13.3k
!v13.3k
)
3134
13.2k
      check_gbr_integer_feasible(cgbr);
3135
13.3k
  }
3136
18.1k
  
if (18.1k
update18.1k
)
3137
3.98k
    cgbr->tab = check_samples(cgbr->tab, ineq, 0);
3138
18.1k
  return;
3139
0
error:
3140
0
  isl_tab_free(cgbr->tab);
3141
0
  cgbr->tab = NULL;
3142
0
}
3143
3144
static isl_stat context_gbr_add_ineq_wrap(void *user, isl_int *ineq)
3145
842
{
3146
842
  struct isl_context *context = (struct isl_context *)user;
3147
842
  context_gbr_add_ineq(context, ineq, 0, 0);
3148
842
  return context->op->is_ok(context) ? 
isl_stat_ok842
:
isl_stat_error0
;
3149
842
}
3150
3151
static enum isl_tab_row_sign context_gbr_ineq_sign(struct isl_context *context,
3152
      isl_int *ineq, int strict)
3153
7.82k
{
3154
7.82k
  struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3155
7.82k
  return tab_ineq_sign(cgbr->tab, ineq, strict);
3156
7.82k
}
3157
3158
/* Check whether "ineq" can be added to the tableau without rendering
3159
 * it infeasible.
3160
 */
3161
static int context_gbr_test_ineq(struct isl_context *context, isl_int *ineq)
3162
7.37k
{
3163
7.37k
  struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3164
7.37k
  struct isl_tab_undo *snap;
3165
7.37k
  struct isl_tab_undo *shifted_snap = NULL;
3166
7.37k
  struct isl_tab_undo *cone_snap = NULL;
3167
7.37k
  int feasible;
3168
7.37k
3169
7.37k
  if (!cgbr->tab)
3170
0
    return -1;
3171
7.37k
3172
7.37k
  
if (7.37k
isl_tab_extend_cons(cgbr->tab, 1) < 07.37k
)
3173
0
    return -1;
3174
7.37k
3175
7.37k
  snap = isl_tab_snap(cgbr->tab);
3176
7.37k
  if (cgbr->shifted)
3177
16
    shifted_snap = isl_tab_snap(cgbr->shifted);
3178
7.37k
  if (cgbr->cone)
3179
1.79k
    cone_snap = isl_tab_snap(cgbr->cone);
3180
7.37k
  add_gbr_ineq(cgbr, ineq);
3181
7.37k
  check_gbr_integer_feasible(cgbr);
3182
7.37k
  if (!cgbr->tab)
3183
0
    return -1;
3184
7.37k
  feasible = !cgbr->tab->empty;
3185
7.37k
  if (isl_tab_rollback(cgbr->tab, snap) < 0)
3186
0
    return -1;
3187
7.37k
  
if (7.37k
shifted_snap7.37k
)
{16
3188
16
    if (isl_tab_rollback(cgbr->shifted, shifted_snap))
3189
0
      return -1;
3190
7.35k
  } else 
if (7.35k
cgbr->shifted7.35k
)
{46
3191
46
    isl_tab_free(cgbr->shifted);
3192
46
    cgbr->shifted = NULL;
3193
46
  }
3194
7.37k
  
if (7.37k
cone_snap7.37k
)
{1.79k
3195
1.79k
    if (isl_tab_rollback(cgbr->cone, cone_snap))
3196
0
      return -1;
3197
5.57k
  } else 
if (5.57k
cgbr->cone5.57k
)
{239
3198
239
    isl_tab_free(cgbr->cone);
3199
239
    cgbr->cone = NULL;
3200
239
  }
3201
7.37k
3202
7.37k
  return feasible;
3203
7.37k
}
3204
3205
/* Return the column of the last of the variables associated to
3206
 * a column that has a non-zero coefficient.
3207
 * This function is called in a context where only coefficients
3208
 * of parameters or divs can be non-zero.
3209
 */
3210
static int last_non_zero_var_col(struct isl_tab *tab, isl_int *p)
3211
372
{
3212
372
  int i;
3213
372
  int col;
3214
372
3215
372
  if (tab->n_var == 0)
3216
0
    return -1;
3217
372
3218
804
  
for (i = tab->n_var - 1; 372
i >= 0804
;
--i432
)
{804
3219
804
    if (
i >= tab->n_param && 804
i < tab->n_var - tab->n_div791
)
3220
24
      continue;
3221
780
    
if (780
tab->var[i].is_row780
)
3222
115
      continue;
3223
665
    col = tab->var[i].index;
3224
665
    if (
!665
isl_int_is_zero665
(p[col]))
3225
372
      return col;
3226
665
  }
3227
372
3228
0
  return -1;
3229
372
}
3230
3231
/* Look through all the recently added equalities in the context
3232
 * to see if we can propagate any of them to the main tableau.
3233
 *
3234
 * The newly added equalities in the context are encoded as pairs
3235
 * of inequalities starting at inequality "first".
3236
 *
3237
 * We tentatively add each of these equalities to the main tableau
3238
 * and if this happens to result in a row with a final coefficient
3239
 * that is one or negative one, we use it to kill a column
3240
 * in the main tableau.  Otherwise, we discard the tentatively
3241
 * added row.
3242
 * This tentative addition of equality constraints turns
3243
 * on the undo facility of the tableau.  Turn it off again
3244
 * at the end, assuming it was turned off to begin with.
3245
 *
3246
 * Return 0 on success and -1 on failure.
3247
 */
3248
static int propagate_equalities(struct isl_context_gbr *cgbr,
3249
  struct isl_tab *tab, unsigned first)
3250
188
{
3251
188
  int i;
3252
188
  struct isl_vec *eq = NULL;
3253
188
  isl_bool needs_undo;
3254
188
3255
188
  needs_undo = isl_tab_need_undo(tab);
3256
188
  if (needs_undo < 0)
3257
0
    goto error;
3258
188
  eq = isl_vec_alloc(tab->mat->ctx, 1 + tab->n_var);
3259
188
  if (!eq)
3260
0
    goto error;
3261
188
3262
188
  
if (188
isl_tab_extend_cons(tab, (cgbr->tab->bmap->n_ineq - first)/2) < 0188
)
3263
0
    goto error;
3264
188
3265
188
  isl_seq_clr(eq->el + 1 + tab->n_param,
3266
188
        tab->n_var - tab->n_param - tab->n_div);
3267
560
  for (i = first; 
i < cgbr->tab->bmap->n_ineq560
;
i += 2372
)
{372
3268
372
    int j;
3269
372
    int r;
3270
372
    struct isl_tab_undo *snap;
3271
372
    snap = isl_tab_snap(tab);
3272
372
3273
372
    isl_seq_cpy(eq->el, cgbr->tab->bmap->ineq[i], 1 + tab->n_param);
3274
372
    isl_seq_cpy(eq->el + 1 + tab->n_var - tab->n_div,
3275
372
          cgbr->tab->bmap->ineq[i] + 1 + tab->n_param,
3276
372
          tab->n_div);
3277
372
3278
372
    r = isl_tab_add_row(tab, eq->el);
3279
372
    if (r < 0)
3280
0
      goto error;
3281
372
    r = tab->con[r].index;
3282
372
    j = last_non_zero_var_col(tab, tab->mat->row[r] + 2 + tab->M);
3283
372
    if (
j < 0 || 372
j < tab->n_dead372
||
3284
372
        
!372
isl_int_is_one372
(tab->mat->row[r][0]) ||
3285
372
        
(!372
isl_int_is_one372
(tab->mat->row[r][2 + tab->M + j]) &&
3286
322
         
!322
isl_int_is_negone322
(tab->mat->row[r][2 + tab->M + j])))
{268
3287
268
      if (isl_tab_rollback(tab, snap) < 0)
3288
0
        goto error;
3289
268
      continue;
3290
268
    }
3291
104
    
if (104
isl_tab_pivot(tab, r, j) < 0104
)
3292
0
      goto error;
3293
104
    
if (104
isl_tab_kill_col(tab, j) < 0104
)
3294
0
      goto error;
3295
104
3296
104
    
if (104
restore_lexmin(tab) < 0104
)
3297
0
      goto error;
3298
104
  }
3299
188
3300
188
  
if (188
!needs_undo188
)
3301
188
    isl_tab_clear_undo(tab);
3302
188
  isl_vec_free(eq);
3303
188
3304
188
  return 0;
3305
0
error:
3306
0
  isl_vec_free(eq);
3307
0
  isl_tab_free(cgbr->tab);
3308
0
  cgbr->tab = NULL;
3309
0
  return -1;
3310
188
}
3311
3312
static int context_gbr_detect_equalities(struct isl_context *context,
3313
  struct isl_tab *tab)
3314
370
{
3315
370
  struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3316
370
  unsigned n_ineq;
3317
370
3318
370
  if (
!cgbr->cone370
)
{102
3319
102
    struct isl_basic_set *bset = isl_tab_peek_bset(cgbr->tab);
3320
102
    cgbr->cone = isl_tab_from_recession_cone(bset, 0);
3321
102
    if (!cgbr->cone)
3322
0
      goto error;
3323
102
    
if (102
isl_tab_track_bset(cgbr->cone,102
3324
102
          isl_basic_set_copy(bset)) < 0)
3325
0
      goto error;
3326
102
  }
3327
370
  
if (370
isl_tab_detect_implicit_equalities(cgbr->cone) < 0370
)
3328
0
    goto error;
3329
370
3330
370
  n_ineq = cgbr->tab->bmap->n_ineq;
3331
370
  cgbr->tab = isl_tab_detect_equalities(cgbr->tab, cgbr->cone);
3332
370
  if (!cgbr->tab)
3333
0
    return -1;
3334
370
  
if (370
cgbr->tab->bmap->n_ineq > n_ineq &&370
3335
188
      propagate_equalities(cgbr, tab, n_ineq) < 0)
3336
0
    return -1;
3337
370
3338
370
  return 0;
3339
0
error:
3340
0
  isl_tab_free(cgbr->tab);
3341
0
  cgbr->tab = NULL;
3342
0
  return -1;
3343
370
}
3344
3345
static int context_gbr_get_div(struct isl_context *context, struct isl_tab *tab,
3346
    struct isl_vec *div)
3347
448
{
3348
448
  return get_div(tab, context, div);
3349
448
}
3350
3351
static isl_bool context_gbr_insert_div(struct isl_context *context, int pos,
3352
  __isl_keep isl_vec *div)
3353
421
{
3354
421
  struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3355
421
  if (
cgbr->cone421
)
{275
3356
275
    int r, n_div, o_div;
3357
275
3358
275
    n_div = isl_basic_map_dim(cgbr->cone->bmap, isl_dim_div);
3359
275
    o_div = cgbr->cone->n_var - n_div;
3360
275
3361
275
    if (isl_tab_extend_cons(cgbr->cone, 3) < 0)
3362
0
      return isl_bool_error;
3363
275
    
if (275
isl_tab_extend_vars(cgbr->cone, 1) < 0275
)
3364
0
      return isl_bool_error;
3365
275
    
if (275
(r = isl_tab_insert_var(cgbr->cone, pos)) <0275
)
3366
0
      return isl_bool_error;
3367
275
3368
275
    cgbr->cone->bmap = isl_basic_map_insert_div(cgbr->cone->bmap,
3369
275
                r - o_div, div);
3370
275
    if (!cgbr->cone->bmap)
3371
0
      return isl_bool_error;
3372
275
    
if (275
isl_tab_push_var(cgbr->cone, isl_tab_undo_bmap_div,275
3373
275
            &cgbr->cone->var[r]) < 0)
3374
0
      return isl_bool_error;
3375
275
  }
3376
421
  return context_tab_insert_div(cgbr->tab, pos, div,
3377
421
          context_gbr_add_ineq_wrap, context);
3378
421
}
3379
3380
static int context_gbr_best_split(struct isl_context *context,
3381
    struct isl_tab *tab)
3382
303
{
3383
303
  struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3384
303
  struct isl_tab_undo *snap;
3385
303
  int r;
3386
303
3387
303
  snap = isl_tab_snap(cgbr->tab);
3388
303
  r = best_split(tab, cgbr->tab);
3389
303
3390
303
  if (
r >= 0 && 303
isl_tab_rollback(cgbr->tab, snap) < 0303
)
3391
0
    return -1;
3392
303
3393
303
  return r;
3394
303
}
3395
3396
static int context_gbr_is_empty(struct isl_context *context)
3397
30.9k
{
3398
30.9k
  struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3399
30.9k
  if (!cgbr->tab)
3400
0
    return -1;
3401
30.9k
  return cgbr->tab->empty;
3402
30.9k
}
3403
3404
struct isl_gbr_tab_undo {
3405
  struct isl_tab_undo *tab_snap;
3406
  struct isl_tab_undo *shifted_snap;
3407
  struct isl_tab_undo *cone_snap;
3408
};
3409
3410
static void *context_gbr_save(struct isl_context *context)
3411
19.0k
{
3412
19.0k
  struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3413
19.0k
  struct isl_gbr_tab_undo *snap;
3414
19.0k
3415
19.0k
  if (!cgbr->tab)
3416
0
    return NULL;
3417
19.0k
3418
19.0k
  
snap = 19.0k
isl_alloc_type19.0k
(cgbr->tab->mat->ctx, struct isl_gbr_tab_undo);
3419
19.0k
  if (!snap)
3420
0
    return NULL;
3421
19.0k
3422
19.0k
  snap->tab_snap = isl_tab_snap(cgbr->tab);
3423
19.0k
  if (isl_tab_save_samples(cgbr->tab) < 0)
3424
0
    goto error;
3425
19.0k
3426
19.0k
  
if (19.0k
cgbr->shifted19.0k
)
3427
10
    snap->shifted_snap = isl_tab_snap(cgbr->shifted);
3428
19.0k
  else
3429
19.0k
    snap->shifted_snap = NULL;
3430
19.0k
3431
19.0k
  if (cgbr->cone)
3432
723
    snap->cone_snap = isl_tab_snap(cgbr->cone);
3433
19.0k
  else
3434
18.3k
    snap->cone_snap = NULL;
3435
19.0k
3436
19.0k
  return snap;
3437
0
error:
3438
0
  free(snap);
3439
0
  return NULL;
3440
19.0k
}
3441
3442
static void context_gbr_restore(struct isl_context *context, void *save)
3443
16.5k
{
3444
16.5k
  struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3445
16.5k
  struct isl_gbr_tab_undo *snap = (struct isl_gbr_tab_undo *)save;
3446
16.5k
  if (!snap)
3447
0
    goto error;
3448
16.5k
  
if (16.5k
isl_tab_rollback(cgbr->tab, snap->tab_snap) < 016.5k
)
3449
0
    goto error;
3450
16.5k
3451
16.5k
  
if (16.5k
snap->shifted_snap16.5k
)
{7
3452
7
    if (isl_tab_rollback(cgbr->shifted, snap->shifted_snap) < 0)
3453
0
      goto error;
3454
16.5k
  } else 
if (16.5k
cgbr->shifted16.5k
)
{0
3455
0
    isl_tab_free(cgbr->shifted);
3456
0
    cgbr->shifted = NULL;
3457
0
  }
3458
16.5k
3459
16.5k
  
if (16.5k
snap->cone_snap16.5k
)
{710
3460
710
    if (isl_tab_rollback(cgbr->cone, snap->cone_snap) < 0)
3461
0
      goto error;
3462
15.7k
  } else 
if (15.7k
cgbr->cone15.7k
)
{155
3463
155
    isl_tab_free(cgbr->cone);
3464
155
    cgbr->cone = NULL;
3465
155
  }
3466
16.5k
3467
16.5k
  free(snap);
3468
16.5k
3469
16.5k
  return;
3470
0
error:
3471
0
  free(snap);
3472
0
  isl_tab_free(cgbr->tab);
3473
0
  cgbr->tab = NULL;
3474
0
}
3475
3476
static void context_gbr_discard(void *save)
3477
2.54k
{
3478
2.54k
  struct isl_gbr_tab_undo *snap = (struct isl_gbr_tab_undo *)save;
3479
2.54k
  free(snap);
3480
2.54k
}
3481
3482
static int context_gbr_is_ok(struct isl_context *context)
3483
1.33k
{
3484
1.33k
  struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3485
1.33k
  return !!cgbr->tab;
3486
1.33k
}
3487
3488
static void context_gbr_invalidate(struct isl_context *context)
3489
0
{
3490
0
  struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3491
0
  isl_tab_free(cgbr->tab);
3492
0
  cgbr->tab = NULL;
3493
0
}
3494
3495
static __isl_null struct isl_context *context_gbr_free(
3496
  struct isl_context *context)
3497
4.76k
{
3498
4.76k
  struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3499
4.76k
  isl_tab_free(cgbr->tab);
3500
4.76k
  isl_tab_free(cgbr->shifted);
3501
4.76k
  isl_tab_free(cgbr->cone);
3502
4.76k
  free(cgbr);
3503
4.76k
3504
4.76k
  return NULL;
3505
4.76k
}
3506
3507
struct isl_context_op isl_context_gbr_op = {
3508
  context_gbr_detect_nonnegative_parameters,
3509
  context_gbr_peek_basic_set,
3510
  context_gbr_peek_tab,
3511
  context_gbr_add_eq,
3512
  context_gbr_add_ineq,
3513
  context_gbr_ineq_sign,
3514
  context_gbr_test_ineq,
3515
  context_gbr_get_div,
3516
  context_gbr_insert_div,
3517
  context_gbr_detect_equalities,
3518
  context_gbr_best_split,
3519
  context_gbr_is_empty,
3520
  context_gbr_is_ok,
3521
  context_gbr_save,
3522
  context_gbr_restore,
3523
  context_gbr_discard,
3524
  context_gbr_invalidate,
3525
  context_gbr_free,
3526
};
3527
3528
static struct isl_context *isl_context_gbr_alloc(__isl_keep isl_basic_set *dom)
3529
4.76k
{
3530
4.76k
  struct isl_context_gbr *cgbr;
3531
4.76k
3532
4.76k
  if (!dom)
3533
0
    return NULL;
3534
4.76k
3535
4.76k
  
cgbr = 4.76k
isl_calloc_type4.76k
(dom->ctx, struct isl_context_gbr);
3536
4.76k
  if (!cgbr)
3537
0
    return NULL;
3538
4.76k
3539
4.76k
  cgbr->context.op = &isl_context_gbr_op;
3540
4.76k
3541
4.76k
  cgbr->shifted = NULL;
3542
4.76k
  cgbr->cone = NULL;
3543
4.76k
  cgbr->tab = isl_tab_from_basic_set(dom, 1);
3544
4.76k
  cgbr->tab = isl_tab_init_samples(cgbr->tab);
3545
4.76k
  if (!cgbr->tab)
3546
0
    goto error;
3547
4.76k
  check_gbr_integer_feasible(cgbr);
3548
4.76k
3549
4.76k
  return &cgbr->context;
3550
0
error:
3551
0
  cgbr->context.op->free(&cgbr->context);
3552
0
  return NULL;
3553
4.76k
}
3554
3555
/* Allocate a context corresponding to "dom".
3556
 * The representation specific fields are initialized by
3557
 * isl_context_lex_alloc or isl_context_gbr_alloc.
3558
 * The shared "n_unknown" field is initialized to the number
3559
 * of final unknown integer divisions in "dom".
3560
 */
3561
static struct isl_context *isl_context_alloc(__isl_keep isl_basic_set *dom)
3562
4.76k
{
3563
4.76k
  struct isl_context *context;
3564
4.76k
  int first;
3565
4.76k
3566
4.76k
  if (!dom)
3567
0
    return NULL;
3568
4.76k
3569
4.76k
  
if (4.76k
dom->ctx->opt->context == 4.76k
ISL_CONTEXT_LEXMIN4.76k
)
3570
0
    context = isl_context_lex_alloc(dom);
3571
4.76k
  else
3572
4.76k
    context = isl_context_gbr_alloc(dom);
3573
4.76k
3574
4.76k
  if (!context)
3575
0
    return NULL;
3576
4.76k
3577
4.76k
  first = isl_basic_set_first_unknown_div(dom);
3578
4.76k
  if (first < 0)
3579
0
    return context->op->free(context);
3580
4.76k
  context->n_unknown = isl_basic_set_dim(dom, isl_dim_div) - first;
3581
4.76k
3582
4.76k
  return context;
3583
4.76k
}
3584
3585
/* Initialize some common fields of "sol", which keeps track
3586
 * of the solution of an optimization problem on "bmap" over
3587
 * the domain "dom".
3588
 * If "max" is set, then a maximization problem is being solved, rather than
3589
 * a minimization problem, which means that the variables in the
3590
 * tableau have value "M - x" rather than "M + x".
3591
 */
3592
static isl_stat sol_init(struct isl_sol *sol, __isl_keep isl_basic_map *bmap,
3593
  __isl_keep isl_basic_set *dom, int max)
3594
4.76k
{
3595
4.76k
  sol->rational = ISL_F_ISSET(bmap, ISL_BASIC_MAP_RATIONAL);
3596
4.76k
  sol->dec_level.callback.run = &sol_dec_level_wrap;
3597
4.76k
  sol->dec_level.sol = sol;
3598
4.76k
  sol->max = max;
3599
4.76k
  sol->n_out = isl_basic_map_dim(bmap, isl_dim_out);
3600
4.76k
  sol->space = isl_basic_map_get_space(bmap);
3601
4.76k
3602
4.76k
  sol->context = isl_context_alloc(dom);
3603
4.76k
  if (
!sol->space || 4.76k
!sol->context4.76k
)
3604
0
    return isl_stat_error;
3605
4.76k
3606
4.76k
  return isl_stat_ok;
3607
4.76k
}
3608
3609
/* Construct an isl_sol_map structure for accumulating the solution.
3610
 * If track_empty is set, then we also keep track of the parts
3611
 * of the context where there is no solution.
3612
 * If max is set, then we are solving a maximization, rather than
3613
 * a minimization problem, which means that the variables in the
3614
 * tableau have value "M - x" rather than "M + x".
3615
 */
3616
static struct isl_sol *sol_map_init(__isl_keep isl_basic_map *bmap,
3617
  __isl_take isl_basic_set *dom, int track_empty, int max)
3618
2.20k
{
3619
2.20k
  struct isl_sol_map *sol_map = NULL;
3620
2.20k
  isl_space *space;
3621
2.20k
3622
2.20k
  if (!bmap)
3623
0
    goto error;
3624
2.20k
3625
2.20k
  
sol_map = 2.20k
isl_calloc_type2.20k
(bmap->ctx, struct isl_sol_map);
3626
2.20k
  if (!sol_map)
3627
0
    goto error;
3628
2.20k
3629
2.20k
  sol_map->sol.free = &sol_map_free;
3630
2.20k
  if (sol_init(&sol_map->sol, bmap, dom, max) < 0)
3631
0
    goto error;
3632
2.20k
  sol_map->sol.add = &sol_map_add_wrap;
3633
1.86k
  sol_map->sol.add_empty = track_empty ? &sol_map_add_empty_wrap : NULL;
3634
2.20k
  space = isl_space_copy(sol_map->sol.space);
3635
2.20k
  sol_map->map = isl_map_alloc_space(space, 1, ISL_MAP_DISJOINT);
3636
2.20k
  if (!sol_map->map)
3637
0
    goto error;
3638
2.20k
3639
2.20k
  
if (2.20k
track_empty2.20k
)
{1.86k
3640
1.86k
    sol_map->empty = isl_set_alloc_space(isl_basic_set_get_space(dom),
3641
1.86k
              1, ISL_SET_DISJOINT);
3642
1.86k
    if (!sol_map->empty)
3643
0
      goto error;
3644
1.86k
  }
3645
2.20k
3646
2.20k
  isl_basic_set_free(dom);
3647
2.20k
  return &sol_map->sol;
3648
0
error:
3649
0
  isl_basic_set_free(dom);
3650
0
  sol_free(&sol_map->sol);
3651
0
  return NULL;
3652
2.20k
}
3653
3654
/* Check whether all coefficients of (non-parameter) variables
3655
 * are non-positive, meaning that no pivots can be performed on the row.
3656
 */
3657
static int is_critical(struct isl_tab *tab, int row)
3658
7.82k
{
3659
7.82k
  int j;
3660
7.82k
  unsigned off = 2 + tab->M;
3661
7.82k
3662
33.3k
  for (j = tab->n_dead; 
j < tab->n_col33.3k
;
++j25.5k
)
{27.0k
3663
27.0k
    if (tab->col_var[j] >= 0 &&
3664
17.7k
        (tab->col_var[j] < tab->n_param  ||
3665
1.40k
        tab->col_var[j] >= tab->n_var - tab->n_div))
3666
17.7k
      continue;
3667
27.0k
3668
9.34k
    
if (9.34k
isl_int_is_pos9.34k
(tab->mat->row[row][off + j]))
3669
1.50k
      return 0;
3670
9.34k
  }
3671
7.82k
3672
6.31k
  return 1;
3673
7.82k
}
3674
3675
/* Check whether the inequality represented by vec is strict over the integers,
3676
 * i.e., there are no integer values satisfying the constraint with
3677
 * equality.  This happens if the gcd of the coefficients is not a divisor
3678
 * of the constant term.  If so, scale the constraint down by the gcd
3679
 * of the coefficients.
3680
 */
3681
static int is_strict(struct isl_vec *vec)
3682
9.77k
{
3683
9.77k
  isl_int gcd;
3684
9.77k
  int strict = 0;
3685
9.77k
3686
9.77k
  isl_int_init(gcd);
3687
9.77k
  isl_seq_gcd(vec->el + 1, vec->size - 1, &gcd);
3688
9.77k
  if (
!9.77k
isl_int_is_one9.77k
(gcd))
{992
3689
992
    strict = !isl_int_is_divisible_by(vec->el[0], gcd);
3690
992
    isl_int_fdiv_q(vec->el[0], vec->el[0], gcd);
3691
992
    isl_seq_scale_down(vec->el + 1, vec->el + 1, gcd, vec->size-1);
3692
992
  }
3693
9.77k
  isl_int_clear(gcd);
3694
9.77k
3695
9.77k
  return strict;
3696
9.77k
}
3697
3698
/* Determine the sign of the given row of the main tableau.
3699
 * The result is one of
3700
 *  isl_tab_row_pos: always non-negative; no pivot needed
3701
 *  isl_tab_row_neg: always non-positive; pivot
3702
 *  isl_tab_row_any: can be both positive and negative; split
3703
 *
3704
 * We first handle some simple cases
3705
 *  - the row sign may be known already
3706
 *  - the row may be obviously non-negative
3707
 *  - the parametric constant may be equal to that of another row
3708
 *    for which we know the sign.  This sign will be either "pos" or
3709
 *    "any".  If it had been "neg" then we would have pivoted before.
3710
 *
3711
 * If none of these cases hold, we check the value of the row for each
3712
 * of the currently active samples.  Based on the signs of these values
3713
 * we make an initial determination of the sign of the row.
3714
 *
3715
 *  all zero      ->  unk(nown)
3716
 *  all non-negative    ->  pos
3717
 *  all non-positive    ->  neg
3718
 *  both negative and positive  ->  all
3719
 *
3720
 * If we end up with "all", we are done.
3721
 * Otherwise, we perform a check for positive and/or negative
3722
 * values as follows.
3723
 *
3724
 *  samples        neg         unk         pos
3725
 *  <0 ?          Y        N      Y        N
3726
 *              pos    any      pos
3727
 *  >0 ?       Y      N  Y     N
3728
 *        any    neg  any   neg
3729
 *
3730
 * There is no special sign for "zero", because we can usually treat zero
3731
 * as either non-negative or non-positive, whatever works out best.
3732
 * However, if the row is "critical", meaning that pivoting is impossible
3733
 * then we don't want to limp zero with the non-positive case, because
3734
 * then we we would lose the solution for those values of the parameters
3735
 * where the value of the row is zero.  Instead, we treat 0 as non-negative
3736
 * ensuring a split if the row can attain both zero and negative values.
3737
 * The same happens when the original constraint was one that could not
3738
 * be satisfied with equality by any integer values of the parameters.
3739
 * In this case, we normalize the constraint, but then a value of zero
3740
 * for the normalized constraint is actually a positive value for the
3741
 * original constraint, so again we need to treat zero as non-negative.
3742
 * In both these cases, we have the following decision tree instead:
3743
 *
3744
 *  all non-negative    ->  pos
3745
 *  all negative      ->  neg
3746
 *  both negative and non-negative  ->  all
3747
 *
3748
 *  samples        neg                     pos
3749
 *  <0 ?                        Y        N
3750
 *                     any      pos
3751
 *  >=0 ?      Y      N
3752
 *        any    neg
3753
 */
3754
static enum isl_tab_row_sign row_sign(struct isl_tab *tab,
3755
  struct isl_sol *sol, int row)
3756
45.9k
{
3757
45.9k
  struct isl_vec *ineq = NULL;
3758
45.9k
  enum isl_tab_row_sign res = isl_tab_row_unknown;
3759
45.9k
  int critical;
3760
45.9k
  int strict;
3761
45.9k
  int row2;
3762
45.9k
3763
45.9k
  if (tab->row_sign[row] != isl_tab_row_unknown)
3764
25.5k
    return tab->row_sign[row];
3765
20.3k
  
if (20.3k
is_obviously_nonneg(tab, row)20.3k
)
3766
12.4k
    return isl_tab_row_pos;
3767
82.9k
  
for (row2 = tab->n_redundant; 7.86k
row2 < tab->n_row82.9k
;
++row275.0k
)
{75.1k
3768
75.1k
    if (tab->row_sign[row2] == isl_tab_row_unknown)
3769
23.7k
      continue;
3770
51.3k
    
if (51.3k
identical_parameter_line(tab, row, row2)51.3k
)
3771
47
      return tab->row_sign[row2];
3772
51.3k
  }
3773
7.86k
3774
7.82k
  critical = is_critical(tab, row);
3775
7.82k
3776
7.82k
  ineq = get_row_parameter_ineq(tab, row);
3777
7.82k
  if (!ineq)
3778
0
    goto error;
3779
7.82k
3780
7.82k
  strict = is_strict(ineq);
3781
7.82k
3782
7.82k
  res = sol->context->op->ineq_sign(sol->context, ineq->el,
3783
1.50k
            critical || strict);
3784
7.82k
3785
7.82k
  if (
res == isl_tab_row_unknown || 7.82k
res == isl_tab_row_pos7.79k
)
{5.73k
3786
5.73k
    /* test for negative values */
3787
5.73k
    int feasible;
3788
5.73k
    isl_seq_neg(ineq->el, ineq->el, ineq->size);
3789
5.73k
    isl_int_sub_ui(ineq->el[0], ineq->el[0], 1);
3790
5.73k
3791
5.73k
    feasible = sol->context->op->test_ineq(sol->context, ineq->el);
3792
5.73k
    if (feasible < 0)
3793
0
      goto error;
3794
5.73k
    
if (5.73k
!feasible5.73k
)
3795
5.34k
      res = isl_tab_row_pos;
3796
5.73k
    else
3797
389
      
res = (res == isl_tab_row_unknown) ? 389
isl_tab_row_neg4
3798
385
                 : isl_tab_row_any;
3799
5.73k
    if (
res == isl_tab_row_neg5.73k
)
{4
3800
4
      isl_seq_neg(ineq->el, ineq->el, ineq->size);
3801
4
      isl_int_sub_ui(ineq->el[0], ineq->el[0], 1);
3802
4
    }
3803
5.73k
  }
3804
7.82k
3805
7.82k
  
if (7.82k
res == isl_tab_row_neg7.82k
)
{1.63k
3806
1.63k
    /* test for positive values */
3807
1.63k
    int feasible;
3808
1.63k
    if (
!critical && 1.63k
!strict155
)
3809
151
      isl_int_sub_ui(ineq->el[0], ineq->el[0], 1);
3810
1.63k
3811
1.63k
    feasible = sol->context->op->test_ineq(sol->context, ineq->el);
3812
1.63k
    if (feasible < 0)
3813
0
      goto error;
3814
1.63k
    
if (1.63k
feasible1.63k
)
3815
1.54k
      res = isl_tab_row_any;
3816
1.63k
  }
3817
7.82k
3818
7.82k
  isl_vec_free(ineq);
3819
7.82k
  return res;
3820
0
error:
3821
0
  isl_vec_free(ineq);
3822
0
  return isl_tab_row_unknown;
3823
7.82k
}
3824
3825
static void find_solutions(struct isl_sol *sol, struct isl_tab *tab);
3826
3827
/* Find solutions for values of the parameters that satisfy the given
3828
 * inequality.
3829
 *
3830
 * We currently take a snapshot of the context tableau that is reset
3831
 * when we return from this function, while we make a copy of the main
3832
 * tableau, leaving the original main tableau untouched.
3833
 * These are fairly arbitrary choices.  Making a copy also of the context
3834
 * tableau would obviate the need to undo any changes made to it later,
3835
 * while taking a snapshot of the main tableau could reduce memory usage.
3836
 * If we were to switch to taking a snapshot of the main tableau,
3837
 * we would have to keep in mind that we need to save the row signs
3838
 * and that we need to do this before saving the current basis
3839
 * such that the basis has been restore before we restore the row signs.
3840
 */
3841
static void find_in_pos(struct isl_sol *sol, struct isl_tab *tab, isl_int *ineq)
3842
1.95k
{
3843
1.95k
  void *saved;
3844
1.95k
3845
1.95k
  if (!sol->context)
3846
0
    goto error;
3847
1.95k
  saved = sol->context->op->save(sol->context);
3848
1.95k
3849
1.95k
  tab = isl_tab_dup(tab);
3850
1.95k
  if (!tab)
3851
0
    goto error;
3852
1.95k
3853
1.95k
  sol->context->op->add_ineq(sol->context, ineq, 0, 1);
3854
1.95k
3855
1.95k
  find_solutions(sol, tab);
3856
1.95k
3857
1.95k
  if (!sol->error)
3858
1.95k
    sol->context->op->restore(sol->context, saved);
3859
1.95k
  else
3860
0
    sol->context->op->discard(saved);
3861
1.95k
  return;
3862
0
error:
3863
0
  sol->error = 1;
3864
0
}
3865
3866
/* Record the absence of solutions for those values of the parameters
3867
 * that do not satisfy the given inequality with equality.
3868
 */
3869
static void no_sol_in_strict(struct isl_sol *sol,
3870
  struct isl_tab *tab, struct isl_vec *ineq)
3871
13.2k
{
3872
13.2k
  int empty;
3873
13.2k
  void *saved;
3874
13.2k
3875
13.2k
  if (
!sol->context || 13.2k
sol->error13.2k
)
3876
0
    goto error;
3877
13.2k
  saved = sol->context->op->save(sol->context);
3878
13.2k
3879
13.2k
  isl_int_sub_ui(ineq->el[0], ineq->el[0], 1);
3880
13.2k
3881
13.2k
  sol->context->op->add_ineq(sol->context, ineq->el, 1, 0);
3882
13.2k
  if (!sol->context)
3883
0
    goto error;
3884
13.2k
3885
13.2k
  empty = tab->empty;
3886
13.2k
  tab->empty = 1;
3887
13.2k
  sol_add(sol, tab);
3888
13.2k
  tab->empty = empty;
3889
13.2k
3890
13.2k
  isl_int_add_ui(ineq->el[0], ineq->el[0], 1);
3891
13.2k
3892
13.2k
  sol->context->op->restore(sol->context, saved);
3893
13.2k
  return;
3894
0
error:
3895
0
  sol->error = 1;
3896
0
}
3897
3898
/* Reset all row variables that are marked to have a sign that may
3899
 * be both positive and negative to have an unknown sign.
3900
 */
3901
static void reset_any_to_unknown(struct isl_tab *tab)
3902
1.95k
{
3903
1.95k
  int row;
3904
1.95k
3905
16.3k
  for (row = tab->n_redundant; 
row < tab->n_row16.3k
;
++row14.3k
)
{14.3k
3906
14.3k
    if (!isl_tab_var_from_row(tab, row)->is_nonneg)
3907
108
      continue;
3908
14.2k
    
if (14.2k
tab->row_sign[row] == isl_tab_row_any14.2k
)
3909
2.36k
      tab->row_sign[row] = isl_tab_row_unknown;
3910
14.2k
  }
3911
1.95k
}
3912
3913
/* Compute the lexicographic minimum of the set represented by the main
3914
 * tableau "tab" within the context "sol->context_tab".
3915
 * On entry the sample value of the main tableau is lexicographically
3916
 * less than or equal to this lexicographic minimum.
3917
 * Pivots are performed until a feasible point is found, which is then
3918
 * necessarily equal to the minimum, or until the tableau is found to
3919
 * be infeasible.  Some pivots may need to be performed for only some
3920
 * feasible values of the context tableau.  If so, the context tableau
3921
 * is split into a part where the pivot is needed and a part where it is not.
3922
 *
3923
 * Whenever we enter the main loop, the main tableau is such that no
3924
 * "obvious" pivots need to be performed on it, where "obvious" means
3925
 * that the given row can be seen to be negative without looking at
3926
 * the context tableau.  In particular, for non-parametric problems,
3927
 * no pivots need to be performed on the main tableau.
3928
 * The caller of find_solutions is responsible for making this property
3929
 * hold prior to the first iteration of the loop, while restore_lexmin
3930
 * is called before every other iteration.
3931
 *
3932
 * Inside the main loop, we first examine the signs of the rows of
3933
 * the main tableau within the context of the context tableau.
3934
 * If we find a row that is always non-positive for all values of
3935
 * the parameters satisfying the context tableau and negative for at
3936
 * least one value of the parameters, we perform the appropriate pivot
3937
 * and start over.  An exception is the case where no pivot can be
3938
 * performed on the row.  In this case, we require that the sign of
3939
 * the row is negative for all values of the parameters (rather than just
3940
 * non-positive).  This special case is handled inside row_sign, which
3941
 * will say that the row can have any sign if it determines that it can
3942
 * attain both negative and zero values.
3943
 *
3944
 * If we can't find a row that always requires a pivot, but we can find
3945
 * one or more rows that require a pivot for some values of the parameters
3946
 * (i.e., the row can attain both positive and negative signs), then we split
3947
 * the context tableau into two parts, one where we force the sign to be
3948
 * non-negative and one where we force is to be negative.
3949
 * The non-negative part is handled by a recursive call (through find_in_pos).
3950
 * Upon returning from this call, we continue with the negative part and
3951
 * perform the required pivot.
3952
 *
3953
 * If no such rows can be found, all rows are non-negative and we have
3954
 * found a (rational) feasible point.  If we only wanted a rational point
3955
 * then we are done.
3956
 * Otherwise, we check if all values of the sample point of the tableau
3957
 * are integral for the variables.  If so, we have found the minimal
3958
 * integral point and we are done.
3959
 * If the sample point is not integral, then we need to make a distinction
3960
 * based on whether the constant term is non-integral or the coefficients
3961
 * of the parameters.  Furthermore, in order to decide how to handle
3962
 * the non-integrality, we also need to know whether the coefficients
3963
 * of the other columns in the tableau are integral.  This leads
3964
 * to the following table.  The first two rows do not correspond
3965
 * to a non-integral sample point and are only mentioned for completeness.
3966
 *
3967
 *  constant  parameters  other
3968
 *
3969
 *  int   int   int |
3970
 *  int   int   rat | -> no problem
3971
 *
3972
 *  rat   int   int   -> fail
3973
 *
3974
 *  rat   int   rat   -> cut
3975
 *
3976
 *  int   rat   rat |
3977
 *  rat   rat   rat | -> parametric cut
3978
 *
3979
 *  int   rat   int |
3980
 *  rat   rat   int | -> split context
3981
 *
3982
 * If the parametric constant is completely integral, then there is nothing
3983
 * to be done.  If the constant term is non-integral, but all the other
3984
 * coefficient are integral, then there is nothing that can be done
3985
 * and the tableau has no integral solution.
3986
 * If, on the other hand, one or more of the other columns have rational
3987
 * coefficients, but the parameter coefficients are all integral, then
3988
 * we can perform a regular (non-parametric) cut.
3989
 * Finally, if there is any parameter coefficient that is non-integral,
3990
 * then we need to involve the context tableau.  There are two cases here.
3991
 * If at least one other column has a rational coefficient, then we
3992
 * can perform a parametric cut in the main tableau by adding a new
3993
 * integer division in the context tableau.
3994
 * If all other columns have integral coefficients, then we need to
3995
 * enforce that the rational combination of parameters (c + \sum a_i y_i)/m
3996
 * is always integral.  We do this by introducing an integer division
3997
 * q = floor((c + \sum a_i y_i)/m) and stipulating that its argument should
3998
 * always be integral in the context tableau, i.e., m q = c + \sum a_i y_i.
3999
 * Since q is expressed in the tableau as
4000
 *  c + \sum a_i y_i - m q >= 0
4001
 *  -c - \sum a_i y_i + m q + m - 1 >= 0
4002
 * it is sufficient to add the inequality
4003
 *  -c - \sum a_i y_i + m q >= 0
4004
 * In the part of the context where this inequality does not hold, the
4005
 * main tableau is marked as being empty.
4006
 */
4007
static void find_solutions(struct isl_sol *sol, struct isl_tab *tab)
4008
5.76k
{
4009
5.76k
  struct isl_context *context;
4010
5.76k
  int r;
4011
5.76k
4012
5.76k
  if (
!tab || 5.76k
sol->error5.76k
)
4013
0
    goto error;
4014
5.76k
4015
5.76k
  context = sol->context;
4016
5.76k
4017
5.76k
  if (tab->empty)
4018
0
    goto done;
4019
5.76k
  
if (5.76k
context->op->is_empty(context)5.76k
)
4020
0
    goto done;
4021
5.76k
4022
8.59k
  
for (r = 0; 5.76k
r >= 0 && 8.59k
tab8.59k
&&
!tab->empty8.59k
;
r = restore_lexmin(tab)2.83k
)
{6.80k
4023
6.80k
    int flags;
4024
6.80k
    int row;
4025
6.80k
    enum isl_tab_row_sign sgn;
4026
6.80k
    int split = -1;
4027
6.80k
    int n_split = 0;
4028
6.80k
4029
53.8k
    for (row = tab->n_redundant; 
row < tab->n_row53.8k
;
++row47.0k
)
{47.1k
4030
47.1k
      if (!isl_tab_var_from_row(tab, row)->is_nonneg)
4031
1.19k
        continue;
4032
45.9k
      sgn = row_sign(tab, sol, row);
4033
45.9k
      if (!sgn)
4034
0
        goto error;
4035
45.9k
      tab->row_sign[row] = sgn;
4036
45.9k
      if (sgn == isl_tab_row_any)
4037
2.38k
        n_split++;
4038
45.9k
      if (
sgn == isl_tab_row_any && 45.9k
split == -12.38k
)
4039
1.97k
        split = row;
4040
45.9k
      if (sgn == isl_tab_row_neg)
4041
99
        break;
4042
45.9k
    }
4043
6.80k
    
if (6.80k
row < tab->n_row6.80k
)
4044
99
      continue;
4045
6.70k
    
if (6.70k
split != -16.70k
)
{1.95k
4046
1.95k
      struct isl_vec *ineq;
4047
1.95k
      if (n_split != 1)
4048
303
        split = context->op->best_split(context, tab);
4049
1.95k
      if (split < 0)
4050
0
        goto error;
4051
1.95k
      ineq = get_row_parameter_ineq(tab, split);
4052
1.95k
      if (!ineq)
4053
0
        goto error;
4054
1.95k
      is_strict(ineq);
4055
1.95k
      reset_any_to_unknown(tab);
4056
1.95k
      tab->row_sign[split] = isl_tab_row_pos;
4057
1.95k
      sol_inc_level(sol);
4058
1.95k
      find_in_pos(sol, tab, ineq->el);
4059
1.95k
      tab->row_sign[split] = isl_tab_row_neg;
4060
1.95k
      isl_seq_neg(ineq->el, ineq->el, ineq->size);
4061
1.95k
      isl_int_sub_ui(ineq->el[0], ineq->el[0], 1);
4062
1.95k
      if (!sol->error)
4063
1.95k
        context->op->add_ineq(context, ineq->el, 0, 1);
4064
1.95k
      isl_vec_free(ineq);
4065
1.95k
      if (sol->error)
4066
0
        goto error;
4067
1.95k
      continue;
4068
1.95k
    }
4069
4.75k
    
if (4.75k
tab->rational4.75k
)
4070
3
      break;
4071
4.74k
    row = first_non_integer_row(tab, &flags);
4072
4.74k
    if (row < 0)
4073
3.96k
      break;
4074
783
    
if (783
ISL_FL_ISSET783
(flags, I_PAR))
{335
4075
335
      if (
ISL_FL_ISSET335
(flags, I_VAR))
{0
4076
0
        if (isl_tab_mark_empty(tab) < 0)
4077
0
          goto error;
4078
0
        break;
4079
0
      }
4080
335
      row = add_cut(tab, row);
4081
448
    } else 
if (448
ISL_FL_ISSET448
(flags, I_VAR))
{76
4082
76
      struct isl_vec *div;
4083
76
      struct isl_vec *ineq;
4084
76
      int d;
4085
76
      div = get_row_split_div(tab, row);
4086
76
      if (!div)
4087
0
        goto error;
4088
76
      d = context->op->get_div(context, tab, div);
4089
76
      isl_vec_free(div);
4090
76
      if (d < 0)
4091
0
        goto error;
4092
76
      ineq = ineq_for_div(context->op->peek_basic_set(context), d);
4093
76
      if (!ineq)
4094
0
        goto error;
4095
76
      sol_inc_level(sol);
4096
76
      no_sol_in_strict(sol, tab, ineq);
4097
76
      isl_seq_neg(ineq->el, ineq->el, ineq->size);
4098
76
      context->op->add_ineq(context, ineq->el, 1, 1);
4099
76
      isl_vec_free(ineq);
4100
76
      if (
sol->error || 76
!context->op->is_ok(context)76
)
4101
0
        goto error;
4102
76
      tab = set_row_cst_to_div(tab, row, d);
4103
76
      if (context->op->is_empty(context))
4104
0
        break;
4105
76
    } else
4106
372
      row = add_parametric_cut(tab, row, context);
4107
783
    
if (783
row < 0783
)
4108
0
      goto error;
4109
783
  }
4110
5.76k
  
if (5.76k
r < 05.76k
)
4111
0
    goto error;
4112
5.76k
done:
4113
5.76k
  sol_add(sol, tab);
4114
5.76k
  isl_tab_free(tab);
4115
5.76k
  return;
4116
0
error:
4117
0
  isl_tab_free(tab);
4118
0
  sol->error = 1;
4119
0
}
4120
4121
/* Does "sol" contain a pair of partial solutions that could potentially
4122
 * be merged?
4123
 *
4124
 * We currently only check that "sol" is not in an error state
4125
 * and that there are at least two partial solutions of which the final two
4126
 * are defined at the same level.
4127
 */
4128
static int sol_has_mergeable_solutions(struct isl_sol *sol)
4129
3.80k
{
4130
3.80k
  if (sol->error)
4131
0
    return 0;
4132
3.80k
  
if (3.80k
!sol->partial3.80k
)
4133
59
    return 0;
4134
3.74k
  
if (3.74k
!sol->partial->next3.74k
)
4135
2.41k
    return 0;
4136
1.33k
  return sol->partial->level == sol->partial->next->level;
4137
3.74k
}
4138
4139
/* Compute the lexicographic minimum of the set represented by the main
4140
 * tableau "tab" within the context "sol->context_tab".
4141
 *
4142
 * As a preprocessing step, we first transfer all the purely parametric
4143
 * equalities from the main tableau to the context tableau, i.e.,
4144
 * parameters that have been pivoted to a row.
4145
 * These equalities are ignored by the main algorithm, because the
4146
 * corresponding rows may not be marked as being non-negative.
4147
 * In parts of the context where the added equality does not hold,
4148
 * the main tableau is marked as being empty.
4149
 *
4150
 * Before we embark on the actual computation, we save a copy
4151
 * of the context.  When we return, we check if there are any
4152
 * partial solutions that can potentially be merged.  If so,
4153
 * we perform a rollback to the initial state of the context.
4154
 * The merging of partial solutions happens inside calls to
4155
 * sol_dec_level that are pushed onto the undo stack of the context.
4156
 * If there are no partial solutions that can potentially be merged
4157
 * then the rollback is skipped as it would just be wasted effort.
4158
 */
4159
static void find_solutions_main(struct isl_sol *sol, struct isl_tab *tab)
4160
3.80k
{
4161
3.80k
  int row;
4162
3.80k
  void *saved;
4163
3.80k
4164
3.80k
  if (!tab)
4165
0
    goto error;
4166
3.80k
4167
3.80k
  sol->level = 0;
4168
3.80k
4169
49.3k
  for (row = tab->n_redundant; 
row < tab->n_row49.3k
;
++row45.5k
)
{45.5k
4170
45.5k
    int p;
4171
45.5k
    struct isl_vec *eq;
4172
45.5k
4173
45.5k
    if (tab->row_var[row] < 0)
4174
7.70k
      continue;
4175
37.8k
    
if (37.8k
tab->row_var[row] >= tab->n_param &&37.8k
4176
31.2k
        tab->row_var[row] < tab->n_var - tab->n_div)
4177
31.2k
      continue;
4178
6.60k
    
if (6.60k
tab->row_var[row] < tab->n_param6.60k
)
4179
6.60k
      p = tab->row_var[row];
4180
6.60k
    else
4181
0
      p = tab->row_var[row]
4182
0
        + tab->n_param - (tab->n_var - tab->n_div);
4183
6.60k
4184
6.60k
    eq = isl_vec_alloc(tab->mat->ctx, 1+tab->n_param+tab->n_div);
4185
6.60k
    if (!eq)
4186
0
      goto error;
4187
6.60k
    get_row_parameter_line(tab, row, eq->el);
4188
6.60k
    isl_int_neg(eq->el[1 + p], tab->mat->row[row][0]);
4189
6.60k
    eq = isl_vec_normalize(eq);
4190
6.60k
4191
6.60k
    sol_inc_level(sol);
4192
6.60k
    no_sol_in_strict(sol, tab, eq);
4193
6.60k
4194
6.60k
    isl_seq_neg(eq->el, eq->el, eq->size);
4195
6.60k
    sol_inc_level(sol);
4196
6.60k
    no_sol_in_strict(sol, tab, eq);
4197
6.60k
    isl_seq_neg(eq->el, eq->el, eq->size);
4198
6.60k
4199
6.60k
    sol->context->op->add_eq(sol->context, eq->el, 1, 1);
4200
6.60k
4201
6.60k
    isl_vec_free(eq);
4202
6.60k
4203
6.60k
    if (isl_tab_mark_redundant(tab, row) < 0)
4204
0
      goto error;
4205
6.60k
4206
6.60k
    
if (6.60k
sol->context->op->is_empty(sol->context)6.60k
)
4207
0
      break;
4208
6.60k
4209
6.60k
    row = tab->n_redundant - 1;
4210
6.60k
  }
4211
3.80k
4212
3.80k
  saved = sol->context->op->save(sol->context);
4213
3.80k
4214
3.80k
  find_solutions(sol, tab);
4215
3.80k
4216
3.80k
  if (sol_has_mergeable_solutions(sol))
4217
1.26k
    sol->context->op->restore(sol->context, saved);
4218
3.80k
  else
4219
2.54k
    sol->context->op->discard(saved);
4220
3.80k
4221
3.80k
  sol->level = 0;
4222
3.80k
  sol_pop(sol);
4223
3.80k
4224
3.80k
  return;
4225
0
error:
4226
0
  isl_tab_free(tab);
4227
0
  sol->error = 1;
4228
0
}
4229
4230
/* Check if integer division "div" of "dom" also occurs in "bmap".
4231
 * If so, return its position within the divs.
4232
 * If not, return -1.
4233
 */
4234
static int find_context_div(struct isl_basic_map *bmap,
4235
  struct isl_basic_set *dom, unsigned div)
4236
500
{
4237
500
  int i;
4238
500
  unsigned b_dim = isl_space_dim(bmap->dim, isl_dim_all);
4239
500
  unsigned d_dim = isl_space_dim(dom->dim, isl_dim_all);
4240
500
4241
500
  if (isl_int_is_zero(dom->div[div][0]))
4242
10
    return -1;
4243
490
  
if (490
isl_seq_first_non_zero(dom->div[div] + 2 + d_dim, dom->n_div) != -1490
)
4244
2
    return -1;
4245
490
4246
576
  
for (i = 0; 488
i < bmap->n_div576
;
++i88
)
{144
4247
144
    if (isl_int_is_zero(bmap->div[i][0]))
4248
12
      continue;
4249
132
    
if (132
isl_seq_first_non_zero(bmap->div[i] + 2 + d_dim,132
4250
132
             (b_dim - d_dim) + bmap->n_div) != -1)
4251
18
      continue;
4252
114
    
if (114
isl_seq_eq(bmap->div[i], dom->div[div], 2 + d_dim)114
)
4253
56
      return i;
4254
114
  }
4255
432
  return -1;
4256
488
}
4257
4258
/* The correspondence between the variables in the main tableau,
4259
 * the context tableau, and the input map and domain is as follows.
4260
 * The first n_param and the last n_div variables of the main tableau
4261
 * form the variables of the context tableau.
4262
 * In the basic map, these n_param variables correspond to the
4263
 * parameters and the input dimensions.  In the domain, they correspond
4264
 * to the parameters and the set dimensions.
4265
 * The n_div variables correspond to the integer divisions in the domain.
4266
 * To ensure that everything lines up, we may need to copy some of the
4267
 * integer divisions of the domain to the map.  These have to be placed
4268
 * in the same order as those in the context and they have to be placed
4269
 * after any other integer divisions that the map may have.
4270
 * This function performs the required reordering.
4271
 */
4272
static __isl_give isl_basic_map *align_context_divs(
4273
  __isl_take isl_basic_map *bmap, __isl_keep isl_basic_set *dom)
4274
241
{
4275
241
  int i;
4276
241
  int common = 0;
4277
241
  int other;
4278
241
4279
491
  for (i = 0; 
i < dom->n_div491
;
++i250
)
4280
250
    
if (250
find_context_div(bmap, dom, i) != -1250
)
4281
28
      common++;
4282
241
  other = bmap->n_div - common;
4283
241
  if (
dom->n_div - common > 0241
)
{217
4284
217
    bmap = isl_basic_map_extend_space(bmap, isl_space_copy(bmap->dim),
4285
217
        dom->n_div - common, 0, 0);
4286
217
    if (!bmap)
4287
0
      return NULL;
4288
217
  }
4289
491
  
for (i = 0; 241
i < dom->n_div491
;
++i250
)
{250
4290
250
    int pos = find_context_div(bmap, dom, i);
4291
250
    if (
pos < 0250
)
{222
4292
222
      pos = isl_basic_map_alloc_div(bmap);
4293
222
      if (pos < 0)
4294
0
        goto error;
4295
222
      
isl_int_set_si222
(bmap->div[pos][0], 0);222
4296
222
    }
4297
250
    
if (250
pos != other + i250
)
4298
15
      isl_basic_map_swap_div(bmap, pos, other + i);
4299
250
  }
4300
241
  return bmap;
4301
0
error:
4302
0
  isl_basic_map_free(bmap);
4303
0
  return NULL;
4304
241
}
4305
4306
/* Base case of isl_tab_basic_map_partial_lexopt, after removing
4307
 * some obvious symmetries.
4308
 *
4309
 * We make sure the divs in the domain are properly ordered,
4310
 * because they will be added one by one in the given order
4311
 * during the construction of the solution map.
4312
 * Furthermore, make sure that the known integer divisions
4313
 * appear before any unknown integer division because the solution
4314
 * may depend on the known integer divisions, while anything that
4315
 * depends on any variable starting from the first unknown integer
4316
 * division is ignored in sol_pma_add.
4317
 */
4318
static struct isl_sol *basic_map_partial_lexopt_base_sol(
4319
  __isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom,
4320
  __isl_give isl_set **empty, int max,
4321
  struct isl_sol *(*init)(__isl_keep isl_basic_map *bmap,
4322
        __isl_take isl_basic_set *dom, int track_empty, int max))
4323
4.58k
{
4324
4.58k
  struct isl_tab *tab;
4325
4.58k
  struct isl_sol *sol = NULL;
4326
4.58k
  struct isl_context *context;
4327
4.58k
4328
4.58k
  if (
dom->n_div4.58k
)
{241
4329
241
    dom = isl_basic_set_sort_divs(dom);
4330
241
    bmap = align_context_divs(bmap, dom);
4331
241
  }
4332
4.58k
  sol = init(bmap, dom, !!empty, max);
4333
4.58k
  if (!sol)
4334
0
    goto error;
4335
4.58k
4336
4.58k
  context = sol->context;
4337
4.58k
  if (isl_basic_set_plain_is_empty(context->op->peek_basic_set(context)))
4338
0
    /* nothing */;
4339
4.58k
  else 
if (4.58k
isl_basic_map_plain_is_empty(bmap)4.58k
)
{960
4340
960
    if (sol->add_empty)
4341
960
      sol->add_empty(sol,
4342
960
        isl_basic_set_copy(context->op->peek_basic_set(context)));
4343
3.62k
  } else {
4344
3.62k
    tab = tab_for_lexmin(bmap,
4345
3.62k
            context->op->peek_basic_set(context), 1, max);
4346
3.62k
    tab = context->op->detect_nonnegative_parameters(context, tab);
4347
3.62k
    find_solutions_main(sol, tab);
4348
3.62k
  }
4349
4.58k
  if (sol->error)
4350
0
    goto error;
4351
4.58k
4352
4.58k
  isl_basic_map_free(bmap);
4353
4.58k
  return sol;
4354
0
error:
4355
0
  sol_free(sol);
4356
0
  isl_basic_map_free(bmap);
4357
0
  return NULL;
4358
4.58k
}
4359
4360
/* Base case of isl_tab_basic_map_partial_lexopt, after removing
4361
 * some obvious symmetries.
4362
 *
4363
 * We call basic_map_partial_lexopt_base_sol and extract the results.
4364
 */
4365
static __isl_give isl_map *basic_map_partial_lexopt_base(
4366
  __isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom,
4367
  __isl_give isl_set **empty, int max)
4368
2.20k
{
4369
2.20k
  isl_map *result = NULL;
4370
2.20k
  struct isl_sol *sol;
4371
2.20k
  struct isl_sol_map *sol_map;
4372
2.20k
4373
2.20k
  sol = basic_map_partial_lexopt_base_sol(bmap, dom, empty, max,
4374
2.20k
            &sol_map_init);
4375
2.20k
  if (!sol)
4376
0
    return NULL;
4377
2.20k
  sol_map = (struct isl_sol_map *) sol;
4378
2.20k
4379
2.20k
  result = isl_map_copy(sol_map->map);
4380
2.20k
  if (empty)
4381
1.86k
    *empty = isl_set_copy(sol_map->empty);
4382
2.20k
  sol_free(&sol_map->sol);
4383
2.20k
  return result;
4384
2.20k
}
4385
4386
/* Return a count of the number of occurrences of the "n" first
4387
 * variables in the inequality constraints of "bmap".
4388
 */
4389
static __isl_give int *count_occurrences(__isl_keep isl_basic_map *bmap,
4390
  int n)
4391
4.61k
{
4392
4.61k
  int i, j;
4393
4.61k
  isl_ctx *ctx;
4394
4.61k
  int *occurrences;
4395
4.61k
4396
4.61k
  if (!bmap)
4397
0
    return NULL;
4398
4.61k
  ctx = isl_basic_map_get_ctx(bmap);
4399
4.61k
  occurrences = isl_calloc_array(ctx, int, n);
4400
4.61k
  if (!occurrences)
4401
0
    return NULL;
4402
4.61k
4403
14.6k
  
for (i = 0; 4.61k
i < bmap->n_ineq14.6k
;
++i10.0k
)
{10.0k
4404
57.8k
    for (j = 0; 
j < n57.8k
;
++j47.7k
)
{47.7k
4405
47.7k
      if (
!47.7k
isl_int_is_zero47.7k
(bmap->ineq[i][1 + j]))
4406
8.24k
        occurrences[j]++;
4407
47.7k
    }
4408
10.0k
  }
4409
4.61k
4410
4.61k
  return occurrences;
4411
4.61k
}
4412
4413
/* Do all of the "n" variables with non-zero coefficients in "c"
4414
 * occur in exactly a single constraint.
4415
 * "occurrences" is an array of length "n" containing the number
4416
 * of occurrences of each of the variables in the inequality constraints.
4417
 */
4418
static int single_occurrence(int n, isl_int *c, int *occurrences)
4419
4.81k
{
4420
4.81k
  int i;
4421
4.81k
4422
18.1k
  for (i = 0; 
i < n18.1k
;
++i13.3k
)
{14.7k
4423
14.7k
    if (isl_int_is_zero(c[i]))
4424
13.0k
      continue;
4425
1.76k
    
if (1.76k
occurrences[i] != 11.76k
)
4426
1.43k
      return 0;
4427
1.76k
  }
4428
4.81k
4429
3.38k
  return 1;
4430
4.81k
}
4431
4432
/* Do all of the "n" initial variables that occur in inequality constraint
4433
 * "ineq" of "bmap" only occur in that constraint?
4434
 */
4435
static int all_single_occurrence(__isl_keep isl_basic_map *bmap, int ineq,
4436
  int n)
4437
0
{
4438
0
  int i, j;
4439
0
4440
0
  for (i = 0; 
i < n0
;
++i0
)
{0
4441
0
    if (isl_int_is_zero(bmap->ineq[ineq][1 + i]))
4442
0
      continue;
4443
0
    
for (j = 0; 0
j < bmap->n_ineq0
;
++j0
)
{0
4444
0
      if (j == ineq)
4445
0
        continue;
4446
0
      
if (0
!0
isl_int_is_zero0
(bmap->ineq[j][1 + i]))
4447
0
        return 0;
4448
0
    }
4449
0
  }
4450
0
4451
0
  return 1;
4452
0
}
4453
4454
/* Structure used during detection of parallel constraints.
4455
 * n_in: number of "input" variables: isl_dim_param + isl_dim_in
4456
 * n_out: number of "output" variables: isl_dim_out + isl_dim_div
4457
 * val: the coefficients of the output variables
4458
 */
4459
struct isl_constraint_equal_info {
4460
  unsigned n_in;
4461
  unsigned n_out;
4462
  isl_int *val;
4463
};
4464
4465
/* Check whether the coefficients of the output variables
4466
 * of the constraint in "entry" are equal to info->val.
4467
 */
4468
static int constraint_equal(const void *entry, const void *val)
4469
31
{
4470
31
  isl_int **row = (isl_int **)entry;
4471
31
  const struct isl_constraint_equal_info *info = val;
4472
31
4473
31
  return isl_seq_eq((*row) + 1 + info->n_in, info->val, info->n_out);
4474
31
}
4475
4476
/* Check whether "bmap" has a pair of constraints that have
4477
 * the same coefficients for the output variables.
4478
 * Note that the coefficients of the existentially quantified
4479
 * variables need to be zero since the existentially quantified
4480
 * of the result are usually not the same as those of the input.
4481
 * Furthermore, check that each of the input variables that occur
4482
 * in those constraints does not occur in any other constraint.
4483
 * If so, return true and return the row indices of the two constraints
4484
 * in *first and *second.
4485
 */
4486
static isl_bool parallel_constraints(__isl_keep isl_basic_map *bmap,
4487
  int *first, int *second)
4488
4.61k
{
4489
4.61k
  int i;
4490
4.61k
  isl_ctx *ctx;
4491
4.61k
  int *occurrences = NULL;
4492
4.61k
  struct isl_hash_table *table = NULL;
4493
4.61k
  struct isl_hash_table_entry *entry;
4494
4.61k
  struct isl_constraint_equal_info info;
4495
4.61k
  unsigned n_out;
4496
4.61k
  unsigned n_div;
4497
4.61k
4498
4.61k
  ctx = isl_basic_map_get_ctx(bmap);
4499
4.61k
  table = isl_hash_table_alloc(ctx, bmap->n_ineq);
4500
4.61k
  if (!table)
4501
0
    goto error;
4502
4.61k
4503
4.61k
  info.n_in = isl_basic_map_dim(bmap, isl_dim_param) +
4504
4.61k
        isl_basic_map_dim(bmap, isl_dim_in);
4505
4.61k
  occurrences = count_occurrences(bmap, info.n_in);
4506
4.61k
  if (
info.n_in && 4.61k
!occurrences4.14k
)
4507
0
    goto error;
4508
4.61k
  n_out = isl_basic_map_dim(bmap, isl_dim_out);
4509
4.61k
  n_div = isl_basic_map_dim(bmap, isl_dim_div);
4510
4.61k
  info.n_out = n_out + n_div;
4511
14.5k
  for (i = 0; 
i < bmap->n_ineq14.5k
;
++i9.98k
)
{10.0k
4512
10.0k
    uint32_t hash;
4513
10.0k
4514
10.0k
    info.val = bmap->ineq[i] + 1 + info.n_in;
4515
10.0k
    if (isl_seq_first_non_zero(info.val, n_out) < 0)
4516
5.16k
      continue;
4517
4.84k
    
if (4.84k
isl_seq_first_non_zero(info.val + n_out, n_div) >= 04.84k
)
4518
32
      continue;
4519
4.81k
    
if (4.81k
!single_occurrence(info.n_in, bmap->ineq[i] + 1,4.81k
4520
4.81k
          occurrences))
4521
1.43k
      continue;
4522
3.38k
    hash = isl_seq_get_hash(info.val, info.n_out);
4523
3.38k
    entry = isl_hash_table_find(ctx, table, hash,
4524
3.38k
              constraint_equal, &info, 1);
4525
3.38k
    if (!entry)
4526
0
      goto error;
4527
3.38k
    
if (3.38k
entry->data3.38k
)
4528
31
      break;
4529
3.34k
    entry->data = &bmap->ineq[i];
4530
3.34k
  }
4531
4.61k
4532
4.61k
  
if (4.61k
i < bmap->n_ineq4.61k
)
{31
4533
31
    *first = ((isl_int **)entry->data) - bmap->ineq; 
4534
31
    *second = i;
4535
31
  }
4536
4.61k
4537
4.61k
  isl_hash_table_free(ctx, table);
4538
4.61k
  free(occurrences);
4539
4.61k
4540
4.61k
  return i < bmap->n_ineq;
4541
0
error:
4542
0
  isl_hash_table_free(ctx, table);
4543
0
  free(occurrences);
4544
0
  return isl_bool_error;
4545
4.61k
}
4546
4547
/* Given a set of upper bounds in "var", add constraints to "bset"
4548
 * that make the i-th bound smallest.
4549
 *
4550
 * In particular, if there are n bounds b_i, then add the constraints
4551
 *
4552
 *  b_i <= b_j  for j > i
4553
 *  b_i <  b_j  for j < i
4554
 */
4555
static __isl_give isl_basic_set *select_minimum(__isl_take isl_basic_set *bset,
4556
  __isl_keep isl_mat *var, int i)
4557
98
{
4558
98
  isl_ctx *ctx;
4559
98
  int j, k;
4560
98
4561
98
  ctx = isl_mat_get_ctx(var);
4562
98
4563
294
  for (j = 0; 
j < var->n_row294
;
++j196
)
{196
4564
196
    if (j == i)
4565
98
      continue;
4566
98
    k = isl_basic_set_alloc_inequality(bset);
4567
98
    if (k < 0)
4568
0
      goto error;
4569
98
    isl_seq_combine(bset->ineq[k], ctx->one, var->row[j],
4570
98
        ctx->negone, var->row[i], var->n_col);
4571
98
    isl_int_set_si(bset->ineq[k][var->n_col], 0);
4572
98
    if (j < i)
4573
49
      isl_int_sub_ui(bset->ineq[k][0], bset->ineq[k][0], 1);
4574
98
  }
4575
98
4576
98
  bset = isl_basic_set_finalize(bset);
4577
98
4578
98
  return bset;
4579
0
error:
4580
0
  isl_basic_set_free(bset);
4581
0
  return NULL;
4582
98
}
4583
4584
/* Given a set of upper bounds on the last "input" variable m,
4585
 * construct a set that assigns the minimal upper bound to m, i.e.,
4586
 * construct a set that divides the space into cells where one
4587
 * of the upper bounds is smaller than all the others and assign
4588
 * this upper bound to m.
4589
 *
4590
 * In particular, if there are n bounds b_i, then the result
4591
 * consists of n basic sets, each one of the form
4592
 *
4593
 *  m = b_i
4594
 *  b_i <= b_j  for j > i
4595
 *  b_i <  b_j  for j < i
4596
 */
4597
static __isl_give isl_set *set_minimum(__isl_take isl_space *dim,
4598
  __isl_take isl_mat *var)
4599
31
{
4600
31
  int i, k;
4601
31
  isl_basic_set *bset = NULL;
4602
31
  isl_set *set = NULL;
4603
31
4604
31
  if (
!dim || 31
!var31
)
4605
0
    goto error;
4606
31
4607
31
  set = isl_set_alloc_space(isl_space_copy(dim),
4608
31
        var->n_row, ISL_SET_DISJOINT);
4609
31
4610
93
  for (i = 0; 
i < var->n_row93
;
++i62
)
{62
4611
62
    bset = isl_basic_set_alloc_space(isl_space_copy(dim), 0,
4612
62
                 1, var->n_row - 1);
4613
62
    k = isl_basic_set_alloc_equality(bset);
4614
62
    if (k < 0)
4615
0
      goto error;
4616
62
    isl_seq_cpy(bset->eq[k], var->row[i], var->n_col);
4617
62
    isl_int_set_si(bset->eq[k][var->n_col], -1);
4618
62
    bset = select_minimum(bset, var, i);
4619
62
    set = isl_set_add_basic_set(set, bset);
4620
62
  }
4621
31
4622
31
  isl_space_free(dim);
4623
31
  isl_mat_free(var);
4624
31
  return set;
4625
0
error:
4626
0
  isl_basic_set_free(bset);
4627
0
  isl_set_free(set);
4628
0
  isl_space_free(dim);
4629
0
  isl_mat_free(var);
4630
0
  return NULL;
4631
31
}
4632
4633
/* Given that the last input variable of "bmap" represents the minimum
4634
 * of the bounds in "cst", check whether we need to split the domain
4635
 * based on which bound attains the minimum.
4636
 *
4637
 * A split is needed when the minimum appears in an integer division
4638
 * or in an equality.  Otherwise, it is only needed if it appears in
4639
 * an upper bound that is different from the upper bounds on which it
4640
 * is defined.
4641
 */
4642
static isl_bool need_split_basic_map(__isl_keep isl_basic_map *bmap,
4643
  __isl_keep isl_mat *cst)
4644
21
{
4645
21
  int i, j;
4646
21
  unsigned total;
4647
21
  unsigned pos;
4648
21
4649
21
  pos = cst->n_col - 1;
4650
21
  total = isl_basic_map_dim(bmap, isl_dim_all);
4651
21
4652
25
  for (i = 0; 
i < bmap->n_div25
;
++i4
)
4653
8
    
if (8
!8
isl_int_is_zero8
(bmap->div[i][2 + pos]))
4654
4
      return isl_bool_true;
4655
21
4656
45
  
for (i = 0; 17
i < bmap->n_eq45
;
++i28
)
4657
29
    
if (29
!29
isl_int_is_zero29
(bmap->eq[i][1 + pos]))
4658
1
      return isl_bool_true;
4659
17
4660
74
  
for (i = 0; 16
i < bmap->n_ineq74
;
++i58
)
{68
4661
68
    if (isl_int_is_nonneg(bmap->ineq[i][1 + pos]))
4662
34
      continue;
4663
34
    
if (34
!34
isl_int_is_negone34
(bmap->ineq[i][1 + pos]))
4664
0
      return isl_bool_true;
4665
34
    
if (34
isl_seq_first_non_zero(bmap->ineq[i] + 1 + pos + 1,34
4666
34
             total - pos - 1) >= 0)
4667
4
      return isl_bool_true;
4668
34
4669
54
    
for (j = 0; 30
j < cst->n_row54
;
++j24
)
4670
48
      
if (48
isl_seq_eq(bmap->ineq[i], cst->row[j], cst->n_col)48
)
4671
24
        break;
4672
30
    if (j >= cst->n_row)
4673
6
      return isl_bool_true;
4674
30
  }
4675
16
4676
6
  return isl_bool_false;
4677
16
}
4678
4679
/* Given that the last set variable of "bset" represents the minimum
4680
 * of the bounds in "cst", check whether we need to split the domain
4681
 * based on which bound attains the minimum.
4682
 *
4683
 * We simply call need_split_basic_map here.  This is safe because
4684
 * the position of the minimum is computed from "cst" and not
4685
 * from "bmap".
4686
 */
4687
static isl_bool need_split_basic_set(__isl_keep isl_basic_set *bset,
4688
  __isl_keep isl_mat *cst)
4689
6
{
4690
6
  return need_split_basic_map(bset_to_bmap(bset), cst);
4691
6
}
4692
4693
/* Given that the last set variable of "set" represents the minimum
4694
 * of the bounds in "cst", check whether we need to split the domain
4695
 * based on which bound attains the minimum.
4696
 */
4697
static isl_bool need_split_set(__isl_keep isl_set *set, __isl_keep isl_mat *cst)
4698
0
{
4699
0
  int i;
4700
0
4701
0
  for (i = 0; 
i < set->n0
;
++i0
)
{0
4702
0
    isl_bool split;
4703
0
4704
0
    split = need_split_basic_set(set->p[i], cst);
4705
0
    if (
split < 0 || 0
split0
)
4706
0
      return split;
4707
0
  }
4708
0
4709
0
  return isl_bool_false;
4710
0
}
4711
4712
/* Given a set of which the last set variable is the minimum
4713
 * of the bounds in "cst", split each basic set in the set
4714
 * in pieces where one of the bounds is (strictly) smaller than the others.
4715
 * This subdivision is given in "min_expr".
4716
 * The variable is subsequently projected out.
4717
 *
4718
 * We only do the split when it is needed.
4719
 * For example if the last input variable m = min(a,b) and the only
4720
 * constraints in the given basic set are lower bounds on m,
4721
 * i.e., l <= m = min(a,b), then we can simply project out m
4722
 * to obtain l <= a and l <= b, without having to split on whether
4723
 * m is equal to a or b.
4724
 */
4725
static __isl_give isl_set *split(__isl_take isl_set *empty,
4726
  __isl_take isl_set *min_expr, __isl_take isl_mat *cst)
4727
6
{
4728
6
  int n_in;
4729
6
  int i;
4730
6
  isl_space *dim;
4731
6
  isl_set *res;
4732
6
4733
6
  if (
!empty || 6
!min_expr6
||
!cst6
)
4734
0
    goto error;
4735
6
4736
6
  n_in = isl_set_dim(empty, isl_dim_set);
4737
6
  dim = isl_set_get_space(empty);
4738
6
  dim = isl_space_drop_dims(dim, isl_dim_set, n_in - 1, 1);
4739
6
  res = isl_set_empty(dim);
4740
6
4741
12
  for (i = 0; 
i < empty->n12
;
++i6
)
{6
4742
6
    isl_bool split;
4743
6
    isl_set *set;
4744
6
4745
6
    set = isl_set_from_basic_set(isl_basic_set_copy(empty->p[i]));
4746
6
    split = need_split_basic_set(empty->p[i], cst);
4747
6
    if (split < 0)
4748
0
      set = isl_set_free(set);
4749
6
    else 
if (6
split6
)
4750
6
      set = isl_set_intersect(set, isl_set_copy(min_expr));
4751
6
    set = isl_set_remove_dims(set, isl_dim_set, n_in - 1, 1);
4752
6
4753
6
    res = isl_set_union_disjoint(res, set);
4754
6
  }
4755
6
4756
6
  isl_set_free(empty);
4757
6
  isl_set_free(min_expr);
4758
6
  isl_mat_free(cst);
4759
6
  return res;
4760
0
error:
4761
0
  isl_set_free(empty);
4762
0
  isl_set_free(min_expr);
4763
0
  isl_mat_free(cst);
4764
0
  return NULL;
4765
6
}
4766
4767
/* Given a map of which the last input variable is the minimum
4768
 * of the bounds in "cst", split each basic set in the set
4769
 * in pieces where one of the bounds is (strictly) smaller than the others.
4770
 * This subdivision is given in "min_expr".
4771
 * The variable is subsequently projected out.
4772
 *
4773
 * The implementation is essentially the same as that of "split".
4774
 */
4775
static __isl_give isl_map *split_domain(__isl_take isl_map *opt,
4776
  __isl_take isl_set *min_expr, __isl_take isl_mat *cst)
4777
13
{
4778
13
  int n_in;
4779
13
  int i;
4780
13
  isl_space *dim;
4781
13
  isl_map *res;
4782
13
4783
13
  if (
!opt || 13
!min_expr13
||
!cst13
)
4784
0
    goto error;
4785
13
4786
13
  n_in = isl_map_dim(opt, isl_dim_in);
4787
13
  dim = isl_map_get_space(opt);
4788
13
  dim = isl_space_drop_dims(dim, isl_dim_in, n_in - 1, 1);
4789
13
  res = isl_map_empty(dim);
4790
13
4791
28
  for (i = 0; 
i < opt->n28
;
++i15
)
{15
4792
15
    isl_map *map;