Coverage Report

Created: 2018-02-20 23:11

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