Coverage Report

Created: 2017-10-03 07:32

/Users/buildslave/jenkins/sharedspace/clang-stage2-coverage-R@2/llvm/tools/polly/lib/External/isl/isl_tab_pip.c
Line
Count
Source (jump to first uncovered line)
1
/*
2
 * Copyright 2008-2009 Katholieke Universiteit Leuven
3
 * Copyright 2010      INRIA Saclay
4
 * Copyright 2016-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 three implementations of this interface,
177
 * isl_sol_map, which simply collects the solutions in an isl_map
178
 * and (optionally) the parts of the context where there is no solution
179
 * in an isl_set,
180
 * isl_sol_pma, which collects an isl_pw_multi_aff instead, and
181
 * isl_sol_for, which calls a user-defined function for each part of
182
 * the solution.
183
 */
184
struct isl_sol {
185
  int error;
186
  int rational;
187
  int level;
188
  int max;
189
  int n_out;
190
  isl_space *space;
191
  struct isl_context *context;
192
  struct isl_partial_sol *partial;
193
  void (*add)(struct isl_sol *sol,
194
    __isl_take isl_basic_set *dom, __isl_take isl_multi_aff *ma);
195
  void (*add_empty)(struct isl_sol *sol, struct isl_basic_set *bset);
196
  void (*free)(struct isl_sol *sol);
197
  struct isl_sol_callback dec_level;
198
};
199
200
static void sol_free(struct isl_sol *sol)
201
7.57k
{
202
7.57k
  struct isl_partial_sol *partial, *next;
203
7.57k
  if (!sol)
204
0
    return;
205
7.57k
  
for (partial = sol->partial; 7.57k
partial7.57k
;
partial = next0
) {
206
0
    next = partial->next;
207
0
    isl_basic_set_free(partial->dom);
208
0
    isl_multi_aff_free(partial->ma);
209
0
    free(partial);
210
0
  }
211
7.57k
  isl_space_free(sol->space);
212
7.57k
  if (sol->context)
213
7.57k
    sol->context->op->free(sol->context);
214
7.57k
  sol->free(sol);
215
7.57k
  free(sol);
216
7.57k
}
217
218
/* Push a partial solution represented by a domain and function "ma"
219
 * onto the stack of partial solutions.
220
 * If "ma" is NULL, then "dom" represents a part of the domain
221
 * with no solution.
222
 */
223
static void sol_push_sol(struct isl_sol *sol,
224
  __isl_take isl_basic_set *dom, __isl_take isl_multi_aff *ma)
225
9.15k
{
226
9.15k
  struct isl_partial_sol *partial;
227
9.15k
228
9.15k
  if (
sol->error || 9.15k
!dom9.15k
)
229
0
    goto error;
230
9.15k
231
9.15k
  
partial = 9.15k
isl_alloc_type9.15k
(dom->ctx, struct isl_partial_sol);
232
9.15k
  if (!partial)
233
0
    goto error;
234
9.15k
235
9.15k
  partial->level = sol->level;
236
9.15k
  partial->dom = dom;
237
9.15k
  partial->ma = ma;
238
9.15k
  partial->next = sol->partial;
239
9.15k
240
9.15k
  sol->partial = partial;
241
9.15k
242
9.15k
  return;
243
0
error:
244
0
  isl_basic_set_free(dom);
245
0
  isl_multi_aff_free(ma);
246
0
  sol->error = 1;
247
0
}
248
249
/* Check that the final columns of "M", starting at "first", are zero.
250
 */
251
static isl_stat check_final_columns_are_zero(__isl_keep isl_mat *M,
252
  unsigned first)
253
6.75k
{
254
6.75k
  int i;
255
6.75k
  unsigned rows, cols, n;
256
6.75k
257
6.75k
  if (!M)
258
0
    return isl_stat_error;
259
6.75k
  rows = isl_mat_rows(M);
260
6.75k
  cols = isl_mat_cols(M);
261
6.75k
  n = cols - first;
262
28.5k
  for (i = 0; 
i < rows28.5k
;
++i21.7k
)
263
21.7k
    
if (21.7k
isl_seq_first_non_zero(M->row[i] + first, n) != -121.7k
)
264
0
      isl_die(isl_mat_get_ctx(M), isl_error_internal,
265
6.75k
        "final columns should be zero",
266
6.75k
        return isl_stat_error);
267
6.75k
  return isl_stat_ok;
268
6.75k
}
269
270
/* Set the affine expressions in "ma" according to the rows in "M", which
271
 * are defined over the local space "ls".
272
 * The matrix "M" may have extra (zero) columns beyond the number
273
 * of variables in "ls".
274
 */
275
static __isl_give isl_multi_aff *set_from_affine_matrix(
276
  __isl_take isl_multi_aff *ma, __isl_take isl_local_space *ls,
277
  __isl_take isl_mat *M)
278
6.75k
{
279
6.75k
  int i, dim;
280
6.75k
  isl_aff *aff;
281
6.75k
282
6.75k
  if (
!ma || 6.75k
!ls6.75k
||
!M6.75k
)
283
0
    goto error;
284
6.75k
285
6.75k
  dim = isl_local_space_dim(ls, isl_dim_all);
286
6.75k
  if (check_final_columns_are_zero(M, 1 + dim) < 0)
287
0
    goto error;
288
21.7k
  
for (i = 1; 6.75k
i < M->n_row21.7k
;
++i15.0k
) {
289
15.0k
    aff = isl_aff_alloc(isl_local_space_copy(ls));
290
15.0k
    if (
aff15.0k
) {
291
15.0k
      isl_int_set(aff->v->el[0], M->row[0][0]);
292
15.0k
      isl_seq_cpy(aff->v->el + 1, M->row[i], 1 + dim);
293
15.0k
    }
294
15.0k
    aff = isl_aff_normalize(aff);
295
15.0k
    ma = isl_multi_aff_set_aff(ma, i - 1, aff);
296
15.0k
  }
297
6.75k
  isl_local_space_free(ls);
298
6.75k
  isl_mat_free(M);
299
6.75k
300
6.75k
  return ma;
301
0
error:
302
0
  isl_local_space_free(ls);
303
0
  isl_mat_free(M);
304
0
  isl_multi_aff_free(ma);
305
0
  return NULL;
306
6.75k
}
307
308
/* Push a partial solution represented by a domain and mapping M
309
 * onto the stack of partial solutions.
310
 *
311
 * The affine matrix "M" maps the dimensions of the context
312
 * to the output variables.  Convert it into an isl_multi_aff and
313
 * then call sol_push_sol.
314
 *
315
 * Note that the description of the initial context may have involved
316
 * existentially quantified variables, in which case they also appear
317
 * in "dom".  These need to be removed before creating the affine
318
 * expression because an affine expression cannot be defined in terms
319
 * of existentially quantified variables without a known representation.
320
 * Since newly added integer divisions are inserted before these
321
 * existentially quantified variables, they are still in the final
322
 * positions and the corresponding final columns of "M" are zero
323
 * because align_context_divs adds the existentially quantified
324
 * variables of the context to the main tableau without any constraints and
325
 * any equality constraints that are added later on can only serve
326
 * to eliminate these existentially quantified variables.
327
 */
328
static void sol_push_sol_mat(struct isl_sol *sol,
329
  __isl_take isl_basic_set *dom, __isl_take isl_mat *M)
330
6.75k
{
331
6.75k
  isl_local_space *ls;
332
6.75k
  isl_multi_aff *ma;
333
6.75k
  int n_div, n_known;
334
6.75k
335
6.75k
  n_div = isl_basic_set_dim(dom, isl_dim_div);
336
6.75k
  n_known = n_div - sol->context->n_unknown;
337
6.75k
338
6.75k
  ma = isl_multi_aff_alloc(isl_space_copy(sol->space));
339
6.75k
  ls = isl_basic_set_get_local_space(dom);
340
6.75k
  ls = isl_local_space_drop_dims(ls, isl_dim_div,
341
6.75k
          n_known, n_div - n_known);
342
6.75k
  ma = set_from_affine_matrix(ma, ls, M);
343
6.75k
344
6.75k
  if (!ma)
345
0
    dom = isl_basic_set_free(dom);
346
6.75k
  sol_push_sol(sol, dom, ma);
347
6.75k
}
348
349
/* Pop one partial solution from the partial solution stack and
350
 * pass it on to sol->add or sol->add_empty.
351
 */
352
static void sol_pop_one(struct isl_sol *sol)
353
9.05k
{
354
9.05k
  struct isl_partial_sol *partial;
355
9.05k
356
9.05k
  partial = sol->partial;
357
9.05k
  sol->partial = partial->next;
358
9.05k
359
9.05k
  if (partial->ma)
360
6.66k
    sol->add(sol, partial->dom, partial->ma);
361
9.05k
  else
362
2.39k
    sol->add_empty(sol, partial->dom);
363
9.05k
  free(partial);
364
9.05k
}
365
366
/* Return a fresh copy of the domain represented by the context tableau.
367
 */
368
static struct isl_basic_set *sol_domain(struct isl_sol *sol)
369
9.24k
{
370
9.24k
  struct isl_basic_set *bset;
371
9.24k
372
9.24k
  if (sol->error)
373
0
    return NULL;
374
9.24k
375
9.24k
  bset = isl_basic_set_dup(sol->context->op->peek_basic_set(sol->context));
376
9.24k
  bset = isl_basic_set_update_from_tab(bset,
377
9.24k
      sol->context->op->peek_tab(sol->context));
378
9.24k
379
9.24k
  return bset;
380
9.24k
}
381
382
/* Check whether two partial solutions have the same affine expressions.
383
 */
384
static isl_bool same_solution(struct isl_partial_sol *s1,
385
  struct isl_partial_sol *s2)
386
2.07k
{
387
2.07k
  if (!s1->ma != !s2->ma)
388
1.93k
    return isl_bool_false;
389
146
  
if (146
!s1->ma146
)
390
0
    return isl_bool_true;
391
146
392
146
  return isl_multi_aff_plain_is_equal(s1->ma, s2->ma);
393
146
}
394
395
/* Swap the initial two partial solutions in "sol".
396
 *
397
 * That is, go from
398
 *
399
 *  sol->partial = p1; p1->next = p2; p2->next = p3
400
 *
401
 * to
402
 *
403
 *  sol->partial = p2; p2->next = p1; p1->next = p3
404
 */
405
static void swap_initial(struct isl_sol *sol)
406
55
{
407
55
  struct isl_partial_sol *partial;
408
55
409
55
  partial = sol->partial;
410
55
  sol->partial = partial->next;
411
55
  partial->next = partial->next->next;
412
55
  sol->partial->next = partial;
413
55
}
414
415
/* Combine the initial two partial solution of "sol" into
416
 * a partial solution with the current context domain of "sol" and
417
 * the function description of the second partial solution in the list.
418
 * The level of the new partial solution is set to the current level.
419
 *
420
 * That is, the first two partial solutions (D1,M1) and (D2,M2) are
421
 * replaced by (D,M2), where D is the domain of "sol", which is assumed
422
 * to be the union of D1 and D2, while M1 is assumed to be equal to M2
423
 * (at least on D1).
424
 */
425
static isl_stat combine_initial_into_second(struct isl_sol *sol)
426
96
{
427
96
  struct isl_partial_sol *partial;
428
96
  isl_basic_set *bset;
429
96
430
96
  partial = sol->partial;
431
96
432
96
  bset = sol_domain(sol);
433
96
  isl_basic_set_free(partial->next->dom);
434
96
  partial->next->dom = bset;
435
96
  partial->next->level = sol->level;
436
96
437
96
  if (!bset)
438
0
    return isl_stat_error;
439
96
440
96
  sol->partial = partial->next;
441
96
  isl_basic_set_free(partial->dom);
442
96
  isl_multi_aff_free(partial->ma);
443
96
  free(partial);
444
96
445
96
  return isl_stat_ok;
446
96
}
447
448
/* Are "ma1" and "ma2" equal to each other on "dom"?
449
 *
450
 * Combine "ma1" and "ma2" with "dom" and check if the results are the same.
451
 * "dom" may have existentially quantified variables.  Eliminate them first
452
 * as otherwise they would have to be eliminated twice, in a more complicated
453
 * context.
454
 */
455
static isl_bool equal_on_domain(__isl_keep isl_multi_aff *ma1,
456
  __isl_keep isl_multi_aff *ma2, __isl_keep isl_basic_set *dom)
457
227
{
458
227
  isl_set *set;
459
227
  isl_pw_multi_aff *pma1, *pma2;
460
227
  isl_bool equal;
461
227
462
227
  set = isl_basic_set_compute_divs(isl_basic_set_copy(dom));
463
227
  pma1 = isl_pw_multi_aff_alloc(isl_set_copy(set),
464
227
          isl_multi_aff_copy(ma1));
465
227
  pma2 = isl_pw_multi_aff_alloc(set, isl_multi_aff_copy(ma2));
466
227
  equal = isl_pw_multi_aff_is_equal(pma1, pma2);
467
227
  isl_pw_multi_aff_free(pma1);
468
227
  isl_pw_multi_aff_free(pma2);
469
227
470
227
  return equal;
471
227
}
472
473
/* The initial two partial solutions of "sol" are known to be at
474
 * the same level.
475
 * If they represent the same solution (on different parts of the domain),
476
 * then combine them into a single solution at the current level.
477
 * Otherwise, pop them both.
478
 *
479
 * Even if the two partial solution are not obviously the same,
480
 * one may still be a simplification of the other over its own domain.
481
 * Also check if the two sets of affine functions are equal when
482
 * restricted to one of the domains.  If so, combine the two
483
 * using the set of affine functions on the other domain.
484
 * That is, for two partial solutions (D1,M1) and (D2,M2),
485
 * if M1 = M2 on D1, then the pair of partial solutions can
486
 * be replaced by (D1+D2,M2) and similarly when M1 = M2 on D2.
487
 */
488
static isl_stat combine_initial_if_equal(struct isl_sol *sol)
489
2.07k
{
490
2.07k
  struct isl_partial_sol *partial;
491
2.07k
  isl_bool same;
492
2.07k
493
2.07k
  partial = sol->partial;
494
2.07k
495
2.07k
  same = same_solution(partial, partial->next);
496
2.07k
  if (same < 0)
497
0
    return isl_stat_error;
498
2.07k
  
if (2.07k
same2.07k
)
499
24
    return combine_initial_into_second(sol);
500
2.05k
  
if (2.05k
partial->ma && 2.05k
partial->next->ma140
) {
501
122
    same = equal_on_domain(partial->ma, partial->next->ma,
502
122
          partial->dom);
503
122
    if (same < 0)
504
0
      return isl_stat_error;
505
122
    
if (122
same122
)
506
17
      return combine_initial_into_second(sol);
507
105
    same = equal_on_domain(partial->ma, partial->next->ma,
508
105
          partial->next->dom);
509
105
    if (
same105
) {
510
55
      swap_initial(sol);
511
55
      return combine_initial_into_second(sol);
512
55
    }
513
1.98k
  }
514
1.98k
515
1.98k
  sol_pop_one(sol);
516
1.98k
  sol_pop_one(sol);
517
1.98k
518
1.98k
  return isl_stat_ok;
519
1.98k
}
520
521
/* Pop all solutions from the partial solution stack that were pushed onto
522
 * the stack at levels that are deeper than the current level.
523
 * If the two topmost elements on the stack have the same level
524
 * and represent the same solution, then their domains are combined.
525
 * This combined domain is the same as the current context domain
526
 * as sol_pop is called each time we move back to a higher level.
527
 * If the outer level (0) has been reached, then all partial solutions
528
 * at the current level are also popped off.
529
 */
530
static void sol_pop(struct isl_sol *sol)
531
8.93k
{
532
8.93k
  struct isl_partial_sol *partial;
533
8.93k
534
8.93k
  if (sol->error)
535
0
    return;
536
8.93k
537
8.93k
  partial = sol->partial;
538
8.93k
  if (!partial)
539
1.99k
    return;
540
6.93k
541
6.93k
  
if (6.93k
partial->level == 0 && 6.93k
sol->level == 02.93k
) {
542
5.87k
    for (partial = sol->partial; 
partial5.87k
;
partial = sol->partial2.93k
)
543
2.93k
      sol_pop_one(sol);
544
2.93k
    return;
545
2.93k
  }
546
4.00k
547
4.00k
  
if (4.00k
partial->level <= sol->level4.00k
)
548
8
    return;
549
3.99k
550
3.99k
  
if (3.99k
partial->next && 3.99k
partial->next->level == partial->level2.23k
) {
551
2.07k
    if (combine_initial_if_equal(sol) < 0)
552
0
      goto error;
553
3.99k
  } else
554
1.91k
    sol_pop_one(sol);
555
3.99k
556
3.99k
  
if (3.99k
sol->level == 03.99k
) {
557
2.02k
    for (partial = sol->partial; 
partial2.02k
;
partial = sol->partial241
)
558
241
      sol_pop_one(sol);
559
1.78k
    return;
560
1.78k
  }
561
2.20k
562
2.20k
  
if (2.20k
02.20k
)
563
0
error:    sol->error = 1;
564
8.93k
}
565
566
static void sol_dec_level(struct isl_sol *sol)
567
2.40k
{
568
2.40k
  if (sol->error)
569
0
    return;
570
2.40k
571
2.40k
  sol->level--;
572
2.40k
573
2.40k
  sol_pop(sol);
574
2.40k
}
575
576
static isl_stat sol_dec_level_wrap(struct isl_tab_callback *cb)
577
2.40k
{
578
2.40k
  struct isl_sol_callback *callback = (struct isl_sol_callback *)cb;
579
2.40k
580
2.40k
  sol_dec_level(callback->sol);
581
2.40k
582
2.40k
  return callback->sol->error ? 
isl_stat_error0
:
isl_stat_ok2.40k
;
583
2.40k
}
584
585
/* Move down to next level and push callback onto context tableau
586
 * to decrease the level again when it gets rolled back across
587
 * the current state.  That is, dec_level will be called with
588
 * the context tableau in the same state as it is when inc_level
589
 * is called.
590
 */
591
static void sol_inc_level(struct isl_sol *sol)
592
21.1k
{
593
21.1k
  struct isl_tab *tab;
594
21.1k
595
21.1k
  if (sol->error)
596
0
    return;
597
21.1k
598
21.1k
  sol->level++;
599
21.1k
  tab = sol->context->op->peek_tab(sol->context);
600
21.1k
  if (isl_tab_push_callback(tab, &sol->dec_level.callback) < 0)
601
0
    sol->error = 1;
602
21.1k
}
603
604
static void scale_rows(struct isl_mat *mat, isl_int m, int n_row)
605
15.0k
{
606
15.0k
  int i;
607
15.0k
608
15.0k
  if (isl_int_is_one(m))
609
15.0k
    return;
610
18
611
51
  
for (i = 0; 18
i < n_row51
;
++i33
)
612
33
    isl_seq_scale(mat->row[i], mat->row[i], m, mat->n_col);
613
15.0k
}
614
615
/* Add the solution identified by the tableau and the context tableau.
616
 *
617
 * The layout of the variables is as follows.
618
 *  tab->n_var is equal to the total number of variables in the input
619
 *      map (including divs that were copied from the context)
620
 *      + the number of extra divs constructed
621
 *      Of these, the first tab->n_param and the last tab->n_div variables
622
 *  correspond to the variables in the context, i.e.,
623
 *    tab->n_param + tab->n_div = context_tab->n_var
624
 *  tab->n_param is equal to the number of parameters and input
625
 *      dimensions in the input map
626
 *  tab->n_div is equal to the number of divs in the context
627
 *
628
 * If there is no solution, then call add_empty with a basic set
629
 * that corresponds to the context tableau.  (If add_empty is NULL,
630
 * then do nothing).
631
 *
632
 * If there is a solution, then first construct a matrix that maps
633
 * all dimensions of the context to the output variables, i.e.,
634
 * the output dimensions in the input map.
635
 * The divs in the input map (if any) that do not correspond to any
636
 * div in the context do not appear in the solution.
637
 * The algorithm will make sure that they have an integer value,
638
 * but these values themselves are of no interest.
639
 * We have to be careful not to drop or rearrange any divs in the
640
 * context because that would change the meaning of the matrix.
641
 *
642
 * To extract the value of the output variables, it should be noted
643
 * that we always use a big parameter M in the main tableau and so
644
 * the variable stored in this tableau is not an output variable x itself, but
645
 *  x' = M + x (in case of minimization)
646
 * or
647
 *  x' = M - x (in case of maximization)
648
 * If x' appears in a column, then its optimal value is zero,
649
 * which means that the optimal value of x is an unbounded number
650
 * (-M for minimization and M for maximization).
651
 * We currently assume that the output dimensions in the original map
652
 * are bounded, so this cannot occur.
653
 * Similarly, when x' appears in a row, then the coefficient of M in that
654
 * row is necessarily 1.
655
 * If the row in the tableau represents
656
 *  d x' = c + d M + e(y)
657
 * then, in case of minimization, the corresponding row in the matrix
658
 * will be
659
 *  a c + a e(y)
660
 * with a d = m, the (updated) common denominator of the matrix.
661
 * In case of maximization, the row will be
662
 *  -a c - a e(y)
663
 */
664
static void sol_add(struct isl_sol *sol, struct isl_tab *tab)
665
27.7k
{
666
27.7k
  struct isl_basic_set *bset = NULL;
667
27.7k
  struct isl_mat *mat = NULL;
668
27.7k
  unsigned off;
669
27.7k
  int row;
670
27.7k
  isl_int m;
671
27.7k
672
27.7k
  if (
sol->error || 27.7k
!tab27.7k
)
673
0
    goto error;
674
27.7k
675
27.7k
  
if (27.7k
tab->empty && 27.7k
!sol->add_empty20.9k
)
676
2.44k
    return;
677
25.2k
  
if (25.2k
sol->context->op->is_empty(sol->context)25.2k
)
678
16.1k
    return;
679
9.15k
680
9.15k
  bset = sol_domain(sol);
681
9.15k
682
9.15k
  if (
tab->empty9.15k
) {
683
2.39k
    sol_push_sol(sol, bset, NULL);
684
2.39k
    return;
685
2.39k
  }
686
6.75k
687
6.75k
  off = 2 + tab->M;
688
6.75k
689
6.75k
  mat = isl_mat_alloc(tab->mat->ctx, 1 + sol->n_out,
690
6.75k
              1 + tab->n_param + tab->n_div);
691
6.75k
  if (!mat)
692
0
    goto error;
693
6.75k
694
6.75k
  
isl_int_init6.75k
(m);
695
6.75k
696
6.75k
  isl_seq_clr(mat->row[0] + 1, mat->n_col - 1);
697
6.75k
  isl_int_set_si(mat->row[0][0], 1);
698
21.7k
  for (row = 0; 
row < sol->n_out21.7k
;
++row15.0k
) {
699
15.0k
    int i = tab->n_param + row;
700
15.0k
    int r, j;
701
15.0k
702
15.0k
    isl_seq_clr(mat->row[1 + row], mat->n_col);
703
15.0k
    if (
!tab->var[i].is_row15.0k
) {
704
0
      if (tab->M)
705
0
        isl_die(mat->ctx, isl_error_invalid,
706
0
          "unbounded optimum", goto error2);
707
0
      continue;
708
15.0k
    }
709
15.0k
710
15.0k
    r = tab->var[i].index;
711
15.0k
    if (tab->M &&
712
15.0k
        isl_int_ne(tab->mat->row[r][2], tab->mat->row[r][0]))
713
0
      isl_die(mat->ctx, isl_error_invalid,
714
15.0k
        "unbounded optimum", goto error2);
715
15.0k
    
isl_int_gcd15.0k
(m, mat->row[0][0], tab->mat->row[r][0]);
716
15.0k
    isl_int_divexact(m, tab->mat->row[r][0], m);
717
15.0k
    scale_rows(mat, m, 1 + row);
718
15.0k
    isl_int_divexact(m, mat->row[0][0], tab->mat->row[r][0]);
719
15.0k
    isl_int_mul(mat->row[1 + row][0], m, tab->mat->row[r][1]);
720
80.0k
    for (j = 0; 
j < tab->n_param80.0k
;
++j65.0k
) {
721
65.0k
      int col;
722
65.0k
      if (tab->var[j].is_row)
723
35.3k
        continue;
724
29.6k
      col = tab->var[j].index;
725
29.6k
      isl_int_mul(mat->row[1 + row][1 + j], m,
726
65.0k
            tab->mat->row[r][off + col]);
727
65.0k
    }
728
16.8k
    for (j = 0; 
j < tab->n_div16.8k
;
++j1.82k
) {
729
1.82k
      int col;
730
1.82k
      if (tab->var[tab->n_var - tab->n_div+j].is_row)
731
170
        continue;
732
1.65k
      col = tab->var[tab->n_var - tab->n_div+j].index;
733
1.65k
      isl_int_mul(mat->row[1 + row][1 + tab->n_param + j], m,
734
1.82k
            tab->mat->row[r][off + col]);
735
1.82k
    }
736
15.0k
    if (sol->max)
737
11.7k
      isl_seq_neg(mat->row[1 + row], mat->row[1 + row],
738
11.7k
            mat->n_col);
739
15.0k
  }
740
6.75k
741
6.75k
  
isl_int_clear6.75k
(m);
742
6.75k
743
6.75k
  sol_push_sol_mat(sol, bset, mat);
744
6.75k
  return;
745
0
error2:
746
0
  isl_int_clear(m);
747
0
error:
748
0
  isl_basic_set_free(bset);
749
0
  isl_mat_free(mat);
750
0
  sol->error = 1;
751
0
}
752
753
struct isl_sol_map {
754
  struct isl_sol  sol;
755
  struct isl_map  *map;
756
  struct isl_set  *empty;
757
};
758
759
static void sol_map_free(struct isl_sol *sol)
760
3.40k
{
761
3.40k
  struct isl_sol_map *sol_map = (struct isl_sol_map *) sol;
762
3.40k
  isl_map_free(sol_map->map);
763
3.40k
  isl_set_free(sol_map->empty);
764
3.40k
}
765
766
/* This function is called for parts of the context where there is
767
 * no solution, with "bset" corresponding to the context tableau.
768
 * Simply add the basic set to the set "empty".
769
 */
770
static void sol_map_add_empty(struct isl_sol_map *sol,
771
  struct isl_basic_set *bset)
772
2.16k
{
773
2.16k
  if (
!bset || 2.16k
!sol->empty2.16k
)
774
0
    goto error;
775
2.16k
776
2.16k
  sol->empty = isl_set_grow(sol->empty, 1);
777
2.16k
  bset = isl_basic_set_simplify(bset);
778
2.16k
  bset = isl_basic_set_finalize(bset);
779
2.16k
  sol->empty = isl_set_add_basic_set(sol->empty, isl_basic_set_copy(bset));
780
2.16k
  if (!sol->empty)
781
0
    goto error;
782
2.16k
  isl_basic_set_free(bset);
783
2.16k
  return;
784
0
error:
785
0
  isl_basic_set_free(bset);
786
0
  sol->sol.error = 1;
787
0
}
788
789
static void sol_map_add_empty_wrap(struct isl_sol *sol,
790
  struct isl_basic_set *bset)
791
2.16k
{
792
2.16k
  sol_map_add_empty((struct isl_sol_map *)sol, bset);
793
2.16k
}
794
795
/* Given a basic set "dom" that represents the context and a tuple of
796
 * affine expressions "ma" defined over this domain, construct a basic map
797
 * that expresses this function on the domain.
798
 */
799
static void sol_map_add(struct isl_sol_map *sol,
800
  __isl_take isl_basic_set *dom, __isl_take isl_multi_aff *ma)
801
3.22k
{
802
3.22k
  isl_basic_map *bmap;
803
3.22k
804
3.22k
  if (
sol->sol.error || 3.22k
!dom3.22k
||
!ma3.22k
)
805
0
    goto error;
806
3.22k
807
3.22k
  bmap = isl_basic_map_from_multi_aff2(ma, sol->sol.rational);
808
3.22k
  bmap = isl_basic_map_intersect_domain(bmap, dom);
809
3.22k
  sol->map = isl_map_grow(sol->map, 1);
810
3.22k
  sol->map = isl_map_add_basic_map(sol->map, bmap);
811
3.22k
  if (!sol->map)
812
0
    sol->sol.error = 1;
813
3.22k
  return;
814
0
error:
815
0
  isl_basic_set_free(dom);
816
0
  isl_multi_aff_free(ma);
817
0
  sol->sol.error = 1;
818
0
}
819
820
static void sol_map_add_wrap(struct isl_sol *sol,
821
  __isl_take isl_basic_set *dom, __isl_take isl_multi_aff *ma)
822
3.22k
{
823
3.22k
  sol_map_add((struct isl_sol_map *)sol, dom, ma);
824
3.22k
}
825
826
827
/* Store the "parametric constant" of row "row" of tableau "tab" in "line",
828
 * i.e., the constant term and the coefficients of all variables that
829
 * appear in the context tableau.
830
 * Note that the coefficient of the big parameter M is NOT copied.
831
 * The context tableau may not have a big parameter and even when it
832
 * does, it is a different big parameter.
833
 */
834
static void get_row_parameter_line(struct isl_tab *tab, int row, isl_int *line)
835
26.0k
{
836
26.0k
  int i;
837
26.0k
  unsigned off = 2 + tab->M;
838
26.0k
839
26.0k
  isl_int_set(line[0], tab->mat->row[row][1]);
840
152k
  for (i = 0; 
i < tab->n_param152k
;
++i126k
) {
841
126k
    if (tab->var[i].is_row)
842
64.6k
      isl_int_set_si(line[1 + i], 0);
843
62.0k
    else {
844
62.0k
      int col = tab->var[i].index;
845
62.0k
      isl_int_set(line[1 + i], tab->mat->row[row][off + col]);
846
62.0k
    }
847
126k
  }
848
32.2k
  for (i = 0; 
i < tab->n_div32.2k
;
++i6.24k
) {
849
6.24k
    if (tab->var[tab->n_var - tab->n_div + i].is_row)
850
2.02k
      isl_int_set_si(line[1 + tab->n_param + i], 0);
851
4.22k
    else {
852
4.22k
      int col = tab->var[tab->n_var - tab->n_div + i].index;
853
4.22k
      isl_int_set(line[1 + tab->n_param + i],
854
4.22k
            tab->mat->row[row][off + col]);
855
4.22k
    }
856
6.24k
  }
857
26.0k
}
858
859
/* Check if rows "row1" and "row2" have identical "parametric constants",
860
 * as explained above.
861
 * In this case, we also insist that the coefficients of the big parameter
862
 * be the same as the values of the constants will only be the same
863
 * if these coefficients are also the same.
864
 */
865
static int identical_parameter_line(struct isl_tab *tab, int row1, int row2)
866
62.6k
{
867
62.6k
  int i;
868
62.6k
  unsigned off = 2 + tab->M;
869
62.6k
870
62.6k
  if (isl_int_ne(tab->mat->row[row1][1], tab->mat->row[row2][1]))
871
57.9k
    return 0;
872
4.75k
873
4.75k
  
if (4.75k
tab->M && 4.75k
isl_int_ne4.75k
(tab->mat->row[row1][2],
874
4.75k
         tab->mat->row[row2][2]))
875
1.72k
    return 0;
876
3.02k
877
7.52k
  
for (i = 0; 3.02k
i < tab->n_param + tab->n_div7.52k
;
++i4.50k
) {
878
7.37k
    int pos = i < tab->n_param ? i :
879
109
      tab->n_var - tab->n_div + i - tab->n_param;
880
7.48k
    int col;
881
7.48k
882
7.48k
    if (tab->var[pos].is_row)
883
2.83k
      continue;
884
4.64k
    col = tab->var[pos].index;
885
4.64k
    if (isl_int_ne(tab->mat->row[row1][off + col],
886
4.64k
             tab->mat->row[row2][off + col]))
887
2.97k
      return 0;
888
7.48k
  }
889
43
  return 1;
890
62.6k
}
891
892
/* Return an inequality that expresses that the "parametric constant"
893
 * should be non-negative.
894
 * This function is only called when the coefficient of the big parameter
895
 * is equal to zero.
896
 */
897
static struct isl_vec *get_row_parameter_ineq(struct isl_tab *tab, int row)
898
16.3k
{
899
16.3k
  struct isl_vec *ineq;
900
16.3k
901
16.3k
  ineq = isl_vec_alloc(tab->mat->ctx, 1 + tab->n_param + tab->n_div);
902
16.3k
  if (!ineq)
903
0
    return NULL;
904
16.3k
905
16.3k
  get_row_parameter_line(tab, row, ineq->el);
906
16.3k
  if (ineq)
907
16.3k
    ineq = isl_vec_normalize(ineq);
908
16.3k
909
16.3k
  return ineq;
910
16.3k
}
911
912
/* Normalize a div expression of the form
913
 *
914
 *  [(g*f(x) + c)/(g * m)]
915
 *
916
 * with c the constant term and f(x) the remaining coefficients, to
917
 *
918
 *  [(f(x) + [c/g])/m]
919
 */
920
static void normalize_div(__isl_keep isl_vec *div)
921
592
{
922
592
  isl_ctx *ctx = isl_vec_get_ctx(div);
923
592
  int len = div->size - 2;
924
592
925
592
  isl_seq_gcd(div->el + 2, len, &ctx->normalize_gcd);
926
592
  isl_int_gcd(ctx->normalize_gcd, ctx->normalize_gcd, div->el[0]);
927
592
928
592
  if (isl_int_is_one(ctx->normalize_gcd))
929
545
    return;
930
47
931
47
  
isl_int_divexact47
(div->el[0], div->el[0], ctx->normalize_gcd);
932
47
  isl_int_fdiv_q(div->el[1], div->el[1], ctx->normalize_gcd);
933
592
  isl_seq_scale_down(div->el + 2, div->el + 2, ctx->normalize_gcd, len);
934
592
}
935
936
/* Return an integer division for use in a parametric cut based
937
 * on the given row.
938
 * In particular, let the parametric constant of the row be
939
 *
940
 *    \sum_i a_i y_i
941
 *
942
 * where y_0 = 1, but none of the y_i corresponds to the big parameter M.
943
 * The div returned is equal to
944
 *
945
 *    floor(\sum_i {-a_i} y_i) = floor((\sum_i (-a_i mod d) y_i)/d)
946
 */
947
static struct isl_vec *get_row_parameter_div(struct isl_tab *tab, int row)
948
504
{
949
504
  struct isl_vec *div;
950
504
951
504
  div = isl_vec_alloc(tab->mat->ctx, 1 + 1 + tab->n_param + tab->n_div);
952
504
  if (!div)
953
0
    return NULL;
954
504
955
504
  
isl_int_set504
(div->el[0], tab->mat->row[row][0]);
956
504
  get_row_parameter_line(tab, row, div->el + 1);
957
504
  isl_seq_neg(div->el + 1, div->el + 1, div->size - 1);
958
504
  normalize_div(div);
959
504
  isl_seq_fdiv_r(div->el + 1, div->el + 1, div->el[0], div->size - 1);
960
504
961
504
  return div;
962
504
}
963
964
/* Return an integer division for use in transferring an integrality constraint
965
 * to the context.
966
 * In particular, let the parametric constant of the row be
967
 *
968
 *    \sum_i a_i y_i
969
 *
970
 * where y_0 = 1, but none of the y_i corresponds to the big parameter M.
971
 * The the returned div is equal to
972
 *
973
 *    floor(\sum_i {a_i} y_i) = floor((\sum_i (a_i mod d) y_i)/d)
974
 */
975
static struct isl_vec *get_row_split_div(struct isl_tab *tab, int row)
976
88
{
977
88
  struct isl_vec *div;
978
88
979
88
  div = isl_vec_alloc(tab->mat->ctx, 1 + 1 + tab->n_param + tab->n_div);
980
88
  if (!div)
981
0
    return NULL;
982
88
983
88
  
isl_int_set88
(div->el[0], tab->mat->row[row][0]);
984
88
  get_row_parameter_line(tab, row, div->el + 1);
985
88
  normalize_div(div);
986
88
  isl_seq_fdiv_r(div->el + 1, div->el + 1, div->el[0], div->size - 1);
987
88
988
88
  return div;
989
88
}
990
991
/* Construct and return an inequality that expresses an upper bound
992
 * on the given div.
993
 * In particular, if the div is given by
994
 *
995
 *  d = floor(e/m)
996
 *
997
 * then the inequality expresses
998
 *
999
 *  m d <= e
1000
 */
1001
static __isl_give isl_vec *ineq_for_div(__isl_keep isl_basic_set *bset,
1002
  unsigned div)
1003
88
{
1004
88
  unsigned total;
1005
88
  unsigned div_pos;
1006
88
  struct isl_vec *ineq;
1007
88
1008
88
  if (!bset)
1009
0
    return NULL;
1010
88
1011
88
  total = isl_basic_set_total_dim(bset);
1012
88
  div_pos = 1 + total - bset->n_div + div;
1013
88
1014
88
  ineq = isl_vec_alloc(bset->ctx, 1 + total);
1015
88
  if (!ineq)
1016
0
    return NULL;
1017
88
1018
88
  isl_seq_cpy(ineq->el, bset->div[div] + 1, 1 + total);
1019
88
  isl_int_neg(ineq->el[div_pos], bset->div[div][0]);
1020
88
  return ineq;
1021
88
}
1022
1023
/* Given a row in the tableau and a div that was created
1024
 * using get_row_split_div and that has been constrained to equality, i.e.,
1025
 *
1026
 *    d = floor(\sum_i {a_i} y_i) = \sum_i {a_i} y_i
1027
 *
1028
 * replace the expression "\sum_i {a_i} y_i" in the row by d,
1029
 * i.e., we subtract "\sum_i {a_i} y_i" and add 1 d.
1030
 * The coefficients of the non-parameters in the tableau have been
1031
 * verified to be integral.  We can therefore simply replace coefficient b
1032
 * by floor(b).  For the coefficients of the parameters we have
1033
 * floor(a_i) = a_i - {a_i}, while for the other coefficients, we have
1034
 * floor(b) = b.
1035
 */
1036
static struct isl_tab *set_row_cst_to_div(struct isl_tab *tab, int row, int div)
1037
88
{
1038
88
  isl_seq_fdiv_q(tab->mat->row[row] + 1, tab->mat->row[row] + 1,
1039
88
      tab->mat->row[row][0], 1 + tab->M + tab->n_col);
1040
88
1041
88
  isl_int_set_si(tab->mat->row[row][0], 1);
1042
88
1043
88
  if (
tab->var[tab->n_var - tab->n_div + div].is_row88
) {
1044
0
    int drow = tab->var[tab->n_var - tab->n_div + div].index;
1045
0
1046
0
    isl_assert(tab->mat->ctx,
1047
0
      isl_int_is_one(tab->mat->row[drow][0]), goto error);
1048
0
    isl_seq_combine(tab->mat->row[row] + 1,
1049
0
      tab->mat->ctx->one, tab->mat->row[row] + 1,
1050
0
      tab->mat->ctx->one, tab->mat->row[drow] + 1,
1051
0
      1 + tab->M + tab->n_col);
1052
88
  } else {
1053
88
    int dcol = tab->var[tab->n_var - tab->n_div + div].index;
1054
88
1055
88
    isl_int_add_ui(tab->mat->row[row][2 + tab->M + dcol],
1056
88
        tab->mat->row[row][2 + tab->M + dcol], 1);
1057
88
  }
1058
88
1059
88
  return tab;
1060
0
error:
1061
0
  isl_tab_free(tab);
1062
0
  return NULL;
1063
88
}
1064
1065
/* Check if the (parametric) constant of the given row is obviously
1066
 * negative, meaning that we don't need to consult the context tableau.
1067
 * If there is a big parameter and its coefficient is non-zero,
1068
 * then this coefficient determines the outcome.
1069
 * Otherwise, we check whether the constant is negative and
1070
 * all non-zero coefficients of parameters are negative and
1071
 * belong to non-negative parameters.
1072
 */
1073
static int is_obviously_neg(struct isl_tab *tab, int row)
1074
189k
{
1075
189k
  int i;
1076
189k
  int col;
1077
189k
  unsigned off = 2 + tab->M;
1078
189k
1079
189k
  if (
tab->M189k
) {
1080
91.9k
    if (isl_int_is_pos(tab->mat->row[row][2]))
1081
47.6k
      return 0;
1082
44.3k
    
if (44.3k
isl_int_is_neg44.3k
(tab->mat->row[row][2]))
1083
0
      return 1;
1084
142k
  }
1085
142k
1086
142k
  
if (142k
isl_int_is_nonneg142k
(tab->mat->row[row][1]))
1087
128k
    return 0;
1088
27.3k
  
for (i = 0; 13.6k
i < tab->n_param27.3k
;
++i13.7k
) {
1089
23.4k
    /* Eliminated parameter */
1090
23.4k
    if (tab->var[i].is_row)
1091
7.43k
      continue;
1092
16.0k
    col = tab->var[i].index;
1093
16.0k
    if (isl_int_is_zero(tab->mat->row[row][off + col]))
1094
6.25k
      continue;
1095
9.74k
    
if (9.74k
!tab->var[i].is_nonneg9.74k
)
1096
8.79k
      return 0;
1097
956
    
if (956
isl_int_is_pos956
(tab->mat->row[row][off + col]))
1098
916
      return 0;
1099
23.4k
  }
1100
4.38k
  
for (i = 0; 3.93k
i < tab->n_div4.38k
;
++i450
) {
1101
619
    if (tab->var[tab->n_var - tab->n_div + i].is_row)
1102
404
      continue;
1103
215
    col = tab->var[tab->n_var - tab->n_div + i].index;
1104
215
    if (isl_int_is_zero(tab->mat->row[row][off + col]))
1105
46
      continue;
1106
169
    
if (169
!tab->var[tab->n_var - tab->n_div + i].is_nonneg169
)
1107
142
      return 0;
1108
27
    
if (27
isl_int_is_pos27
(tab->mat->row[row][off + col]))
1109
27
      return 0;
1110
619
  }
1111
3.76k
  return 1;
1112
189k
}
1113
1114
/* Check if the (parametric) constant of the given row is obviously
1115
 * non-negative, meaning that we don't need to consult the context tableau.
1116
 * If there is a big parameter and its coefficient is non-zero,
1117
 * then this coefficient determines the outcome.
1118
 * Otherwise, we check whether the constant is non-negative and
1119
 * all non-zero coefficients of parameters are positive and
1120
 * belong to non-negative parameters.
1121
 */
1122
static int is_obviously_nonneg(struct isl_tab *tab, int row)
1123
30.0k
{
1124
30.0k
  int i;
1125
30.0k
  int col;
1126
30.0k
  unsigned off = 2 + tab->M;
1127
30.0k
1128
30.0k
  if (
tab->M30.0k
) {
1129
30.0k
    if (isl_int_is_pos(tab->mat->row[row][2]))
1130
12.5k
      return 1;
1131
17.4k
    
if (17.4k
isl_int_is_neg17.4k
(tab->mat->row[row][2]))
1132
0
      return 0;
1133
17.4k
  }
1134
17.4k
1135
17.4k
  
if (17.4k
isl_int_is_neg17.4k
(tab->mat->row[row][1]))
1136
5.76k
    return 0;
1137
44.4k
  
for (i = 0; 11.6k
i < tab->n_param44.4k
;
++i32.7k
) {
1138
38.1k
    /* Eliminated parameter */
1139
38.1k
    if (tab->var[i].is_row)
1140
16.6k
      continue;
1141
21.5k
    col = tab->var[i].index;
1142
21.5k
    if (isl_int_is_zero(tab->mat->row[row][off + col]))
1143
13.5k
      continue;
1144
8.01k
    
if (8.01k
!tab->var[i].is_nonneg8.01k
)
1145
1.75k
      return 0;
1146
6.26k
    
if (6.26k
isl_int_is_neg6.26k
(tab->mat->row[row][off + col]))
1147
3.65k
      return 0;
1148
38.1k
  }
1149
10.6k
  
for (i = 0; 6.27k
i < tab->n_div10.6k
;
++i4.34k
) {
1150
4.52k
    if (tab->var[tab->n_var - tab->n_div + i].is_row)
1151
3.75k
      continue;
1152
773
    col = tab->var[tab->n_var - tab->n_div + i].index;
1153
773
    if (isl_int_is_zero(tab->mat->row[row][off + col]))
1154
585
      continue;
1155
188
    
if (188
!tab->var[tab->n_var - tab->n_div + i].is_nonneg188
)
1156
174
      return 0;
1157
14
    
if (14
isl_int_is_neg14
(tab->mat->row[row][off + col]))
1158
7
      return 0;
1159
4.52k
  }
1160
6.09k
  return 1;
1161
30.0k
}
1162
1163
/* Given a row r and two columns, return the column that would
1164
 * lead to the lexicographically smallest increment in the sample
1165
 * solution when leaving the basis in favor of the row.
1166
 * Pivoting with column c will increment the sample value by a non-negative
1167
 * constant times a_{V,c}/a_{r,c}, with a_{V,c} the elements of column c
1168
 * corresponding to the non-parametric variables.
1169
 * If variable v appears in a column c_v, the a_{v,c} = 1 iff c = c_v,
1170
 * with all other entries in this virtual row equal to zero.
1171
 * If variable v appears in a row, then a_{v,c} is the element in column c
1172
 * of that row.
1173
 *
1174
 * Let v be the first variable with a_{v,c1}/a_{r,c1} != a_{v,c2}/a_{r,c2}.
1175
 * Then if a_{v,c1}/a_{r,c1} < a_{v,c2}/a_{r,c2}, i.e.,
1176
 * a_{v,c2} a_{r,c1} - a_{v,c1} a_{r,c2} > 0, c1 results in the minimal
1177
 * increment.  Otherwise, it's c2.
1178
 */
1179
static int lexmin_col_pair(struct isl_tab *tab,
1180
  int row, int col1, int col2, isl_int tmp)
1181
4.95k
{
1182
4.95k
  int i;
1183
4.95k
  isl_int *tr;
1184
4.95k
1185
4.95k
  tr = tab->mat->row[row] + 2 + tab->M;
1186
4.95k
1187
21.0k
  for (i = tab->n_param; 
i < tab->n_var - tab->n_div21.0k
;
++i16.1k
) {
1188
21.0k
    int s1, s2;
1189
21.0k
    isl_int *r;
1190
21.0k
1191
21.0k
    if (
!tab->var[i].is_row21.0k
) {
1192
7.20k
      if (tab->var[i].index == col1)
1193
701
        return col2;
1194
6.50k
      
if (6.50k
tab->var[i].index == col26.50k
)
1195
464
        return col1;
1196
6.03k
      continue;
1197
6.03k
    }
1198
13.8k
1199
13.8k
    
if (13.8k
tab->var[i].index == row13.8k
)
1200
106
      continue;
1201
13.7k
1202
13.7k
    r = tab->mat->row[tab->var[i].index] + 2 + tab->M;
1203
13.7k
    s1 = isl_int_sgn(r[col1]);
1204
13.7k
    s2 = isl_int_sgn(r[col2]);
1205
13.7k
    if (
s1 == 0 && 13.7k
s2 == 010.4k
)
1206
8.96k
      continue;
1207
4.81k
    
if (4.81k
s1 < s24.81k
)
1208
1.67k
      return col1;
1209
3.14k
    
if (3.14k
s2 < s13.14k
)
1210
881
      return col2;
1211
2.26k
1212
2.26k
    
isl_int_mul2.26k
(tmp, r[col2], tr[col1]);
1213
2.26k
    isl_int_submul(tmp, r[col1], tr[col2]);
1214
2.26k
    if (isl_int_is_pos(tmp))
1215
830
      return col1;
1216
1.43k
    
if (1.43k
isl_int_is_neg1.43k
(tmp))
1217
408
      return col2;
1218
21.0k
  }
1219
0
  return -1;
1220
4.95k
}
1221
1222
/* Given a row in the tableau, find and return the column that would
1223
 * result in the lexicographically smallest, but positive, increment
1224
 * in the sample point.
1225
 * If there is no such column, then return tab->n_col.
1226
 * If anything goes wrong, return -1.
1227
 */
1228
static int lexmin_pivot_col(struct isl_tab *tab, int row)
1229
11.8k
{
1230
11.8k
  int j;
1231
11.8k
  int col = tab->n_col;
1232
11.8k
  isl_int *tr;
1233
11.8k
  isl_int tmp;
1234
11.8k
1235
11.8k
  tr = tab->mat->row[row] + 2 + tab->M;
1236
11.8k
1237
11.8k
  isl_int_init(tmp);
1238
11.8k
1239
88.2k
  for (j = tab->n_dead; 
j < tab->n_col88.2k
;
++j76.3k
) {
1240
76.3k
    if (tab->col_var[j] >= 0 &&
1241
57.4k
        (tab->col_var[j] < tab->n_param  ||
1242
40.8k
        tab->col_var[j] >= tab->n_var - tab->n_div))
1243
18.4k
      continue;
1244
57.9k
1245
57.9k
    
if (57.9k
!57.9k
isl_int_is_pos57.9k
(tr[j]))
1246
44.1k
      continue;
1247
13.7k
1248
13.7k
    
if (13.7k
col == tab->n_col13.7k
)
1249
8.76k
      col = j;
1250
13.7k
    else
1251
4.95k
      col = lexmin_col_pair(tab, row, col, j, tmp);
1252
13.7k
    isl_assert(tab->mat->ctx, col >= 0, goto error);
1253
76.3k
  }
1254
11.8k
1255
11.8k
  
isl_int_clear11.8k
(tmp);
1256
11.8k
  return col;
1257
0
error:
1258
0
  isl_int_clear(tmp);
1259
0
  return -1;
1260
11.8k
}
1261
1262
/* Return the first known violated constraint, i.e., a non-negative
1263
 * constraint that currently has an either obviously negative value
1264
 * or a previously determined to be negative value.
1265
 *
1266
 * If any constraint has a negative coefficient for the big parameter,
1267
 * if any, then we return one of these first.
1268
 */
1269
static int first_neg(struct isl_tab *tab)
1270
40.5k
{
1271
40.5k
  int row;
1272
40.5k
1273
40.5k
  if (tab->M)
1274
246k
    
for (row = tab->n_redundant; 32.3k
row < tab->n_row246k
;
++row213k
) {
1275
218k
      if (!isl_tab_var_from_row(tab, row)->is_nonneg)
1276
53.8k
        continue;
1277
164k
      
if (164k
!164k
isl_int_is_neg164k
(tab->mat->row[row][2]))
1278
160k
        continue;
1279
4.27k
      
if (4.27k
tab->row_sign4.27k
)
1280
4.27k
        tab->row_sign[row] = isl_tab_row_neg;
1281
32.3k
      return row;
1282
32.3k
    }
1283
314k
  
for (row = tab->n_redundant; 36.3k
row < tab->n_row314k
;
++row278k
) {
1284
286k
    if (!isl_tab_var_from_row(tab, row)->is_nonneg)
1285
46.2k
      continue;
1286
239k
    
if (239k
tab->row_sign239k
) {
1287
142k
      if (tab->row_sign[row] == 0 &&
1288
91.9k
          is_obviously_neg(tab, row))
1289
264
        tab->row_sign[row] = isl_tab_row_neg;
1290
142k
      if (tab->row_sign[row] != isl_tab_row_neg)
1291
137k
        continue;
1292
97.8k
    } else 
if (97.8k
!is_obviously_neg(tab, row)97.8k
)
1293
94.3k
      continue;
1294
7.55k
    return row;
1295
7.55k
  }
1296
28.7k
  return -1;
1297
40.5k
}
1298
1299
/* Check whether the invariant that all columns are lexico-positive
1300
 * is satisfied.  This function is not called from the current code
1301
 * but is useful during debugging.
1302
 */
1303
static void check_lexpos(struct isl_tab *tab) __attribute__ ((unused));
1304
static void check_lexpos(struct isl_tab *tab)
1305
0
{
1306
0
  unsigned off = 2 + tab->M;
1307
0
  int col;
1308
0
  int var;
1309
0
  int row;
1310
0
1311
0
  for (col = tab->n_dead; col < tab->n_col; ++col) {
1312
0
    if (tab->col_var[col] >= 0 &&
1313
0
        (tab->col_var[col] < tab->n_param ||
1314
0
         tab->col_var[col] >= tab->n_var - tab->n_div))
1315
0
      continue;
1316
0
    for (var = tab->n_param; var < tab->n_var - tab->n_div; ++var) {
1317
0
      if (!tab->var[var].is_row) {
1318
0
        if (tab->var[var].index == col)
1319
0
          break;
1320
0
        else
1321
0
          continue;
1322
0
      }
1323
0
      row = tab->var[var].index;
1324
0
      if (isl_int_is_zero(tab->mat->row[row][off + col]))
1325
0
        continue;
1326
0
      if (isl_int_is_pos(tab->mat->row[row][off + col]))
1327
0
        break;
1328
0
      fprintf(stderr, "lexneg column %d (row %d)\n",
1329
0
        col, row);
1330
0
    }
1331
0
    if (var >= tab->n_var - tab->n_div)
1332
0
      fprintf(stderr, "zero column %d\n", col);
1333
0
  }
1334
0
}
1335
1336
/* Report to the caller that the given constraint is part of an encountered
1337
 * conflict.
1338
 */
1339
static int report_conflicting_constraint(struct isl_tab *tab, int con)
1340
1.16k
{
1341
1.16k
  return tab->conflict(con, tab->conflict_user);
1342
1.16k
}
1343
1344
/* Given a conflicting row in the tableau, report all constraints
1345
 * involved in the row to the caller.  That is, the row itself
1346
 * (if it represents a constraint) and all constraint columns with
1347
 * non-zero (and therefore negative) coefficients.
1348
 */
1349
static int report_conflict(struct isl_tab *tab, int row)
1350
3.06k
{
1351
3.06k
  int j;
1352
3.06k
  isl_int *tr;
1353
3.06k
1354
3.06k
  if (!tab->conflict)
1355
2.61k
    return 0;
1356
453
1357
453
  
if (453
tab->row_var[row] < 0 &&
1358
453
      report_conflicting_constraint(tab, ~tab->row_var[row]) < 0)
1359
0
    return -1;
1360
453
1361
453
  tr = tab->mat->row[row] + 2 + tab->M;
1362
453
1363
5.31k
  for (j = tab->n_dead; 
j < tab->n_col5.31k
;
++j4.86k
) {
1364
4.86k
    if (tab->col_var[j] >= 0 &&
1365
3.94k
        (tab->col_var[j] < tab->n_param  ||
1366
3.94k
        tab->col_var[j] >= tab->n_var - tab->n_div))
1367
0
      continue;
1368
4.86k
1369
4.86k
    
if (4.86k
!4.86k
isl_int_is_neg4.86k
(tr[j]))
1370
4.12k
      continue;
1371
743
1372
743
    
if (743
tab->col_var[j] < 0 &&
1373
707
        report_conflicting_constraint(tab, ~tab->col_var[j]) < 0)
1374
0
      return -1;
1375
4.86k
  }
1376
453
1377
453
  return 0;
1378
3.06k
}
1379
1380
/* Resolve all known or obviously violated constraints through pivoting.
1381
 * In particular, as long as we can find any violated constraint, we
1382
 * look for a pivoting column that would result in the lexicographically
1383
 * smallest increment in the sample point.  If there is no such column
1384
 * then the tableau is infeasible.
1385
 */
1386
static int restore_lexmin(struct isl_tab *tab) WARN_UNUSED;
1387
static int restore_lexmin(struct isl_tab *tab)
1388
31.8k
{
1389
31.8k
  int row, col;
1390
31.8k
1391
31.8k
  if (!tab)
1392
0
    return -1;
1393
31.8k
  
if (31.8k
tab->empty31.8k
)
1394
26
    return 0;
1395
40.5k
  
while (31.8k
(row = first_neg(tab)) != -140.5k
) {
1396
11.8k
    col = lexmin_pivot_col(tab, row);
1397
11.8k
    if (
col >= tab->n_col11.8k
) {
1398
3.06k
      if (report_conflict(tab, row) < 0)
1399
0
        return -1;
1400
3.06k
      
if (3.06k
isl_tab_mark_empty(tab) < 03.06k
)
1401
0
        return -1;
1402
3.06k
      return 0;
1403
3.06k
    }
1404
8.76k
    
if (8.76k
col < 08.76k
)
1405
0
      return -1;
1406
8.76k
    
if (8.76k
isl_tab_pivot(tab, row, col) < 08.76k
)
1407
0
      return -1;
1408
11.8k
  }
1409
28.7k
  return 0;
1410
31.8k
}
1411
1412
/* Given a row that represents an equality, look for an appropriate
1413
 * pivoting column.
1414
 * In particular, if there are any non-zero coefficients among
1415
 * the non-parameter variables, then we take the last of these
1416
 * variables.  Eliminating this variable in terms of the other
1417
 * variables and/or parameters does not influence the property
1418
 * that all column in the initial tableau are lexicographically
1419
 * positive.  The row corresponding to the eliminated variable
1420
 * will only have non-zero entries below the diagonal of the
1421
 * initial tableau.  That is, we transform
1422
 *
1423
 *    I       I
1424
 *      1   into    a
1425
 *        I         I
1426
 *
1427
 * If there is no such non-parameter variable, then we are dealing with
1428
 * pure parameter equality and we pick any parameter with coefficient 1 or -1
1429
 * for elimination.  This will ensure that the eliminated parameter
1430
 * always has an integer value whenever all the other parameters are integral.
1431
 * If there is no such parameter then we return -1.
1432
 */
1433
static int last_var_col_or_int_par_col(struct isl_tab *tab, int row)
1434
21.7k
{
1435
21.7k
  unsigned off = 2 + tab->M;
1436
21.7k
  int i;
1437
21.7k
1438
77.6k
  for (i = tab->n_var - tab->n_div - 1; 
i >= 0 && 77.6k
i >= tab->n_param76.7k
;
--i55.9k
) {
1439
67.5k
    int col;
1440
67.5k
    if (tab->var[i].is_row)
1441
39.6k
      continue;
1442
27.9k
    col = tab->var[i].index;
1443
27.9k
    if (col <= tab->n_dead)
1444
2.22k
      continue;
1445
25.6k
    
if (25.6k
!25.6k
isl_int_is_zero25.6k
(tab->mat->row[row][off + col]))
1446
11.6k
      return col;
1447
67.5k
  }
1448
26.8k
  
for (i = tab->n_dead; 10.0k
i < tab->n_col26.8k
;
++i16.7k
) {
1449
26.7k
    if (isl_int_is_one(tab->mat->row[row][off + i]))
1450
4.70k
      return i;
1451
22.0k
    
if (22.0k
isl_int_is_negone22.0k
(tab->mat->row[row][off + i]))
1452
5.29k
      return i;
1453
26.7k
  }
1454
83
  return -1;
1455
21.7k
}
1456
1457
/* Add an equality that is known to be valid to the tableau.
1458
 * We first check if we can eliminate a variable or a parameter.
1459
 * If not, we add the equality as two inequalities.
1460
 * In this case, the equality was a pure parameter equality and there
1461
 * is no need to resolve any constraint violations.
1462
 *
1463
 * This function assumes that at least two more rows and at least
1464
 * two more elements in the constraint array are available in the tableau.
1465
 */
1466
static struct isl_tab *add_lexmin_valid_eq(struct isl_tab *tab, isl_int *eq)
1467
21.7k
{
1468
21.7k
  int i;
1469
21.7k
  int r;
1470
21.7k
1471
21.7k
  if (!tab)
1472
0
    return NULL;
1473
21.7k
  r = isl_tab_add_row(tab, eq);
1474
21.7k
  if (r < 0)
1475
0
    goto error;
1476
21.7k
1477
21.7k
  r = tab->con[r].index;
1478
21.7k
  i = last_var_col_or_int_par_col(tab, r);
1479
21.7k
  if (
i < 021.7k
) {
1480
83
    tab->con[r].is_nonneg = 1;
1481
83
    if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r]) < 0)
1482
0
      goto error;
1483
83
    isl_seq_neg(eq, eq, 1 + tab->n_var);
1484
83
    r = isl_tab_add_row(tab, eq);
1485
83
    if (r < 0)
1486
0
      goto error;
1487
83
    tab->con[r].is_nonneg = 1;
1488
83
    if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r]) < 0)
1489
0
      goto error;
1490
21.6k
  } else {
1491
21.6k
    if (isl_tab_pivot(tab, r, i) < 0)
1492
0
      goto error;
1493
21.6k
    
if (21.6k
isl_tab_kill_col(tab, i) < 021.6k
)
1494
0
      goto error;
1495
21.6k
    tab->n_eq++;
1496
21.6k
  }
1497
21.7k
1498
21.7k
  return tab;
1499
0
error:
1500
0
  isl_tab_free(tab);
1501
0
  return NULL;
1502
21.7k
}
1503
1504
/* Check if the given row is a pure constant.
1505
 */
1506
static int is_constant(struct isl_tab *tab, int row)
1507
301
{
1508
301
  unsigned off = 2 + tab->M;
1509
301
1510
301
  return isl_seq_first_non_zero(tab->mat->row[row] + off + tab->n_dead,
1511
301
          tab->n_col - tab->n_dead) == -1;
1512
301
}
1513
1514
/* Add an equality that may or may not be valid to the tableau.
1515
 * If the resulting row is a pure constant, then it must be zero.
1516
 * Otherwise, the resulting tableau is empty.
1517
 *
1518
 * If the row is not a pure constant, then we add two inequalities,
1519
 * each time checking that they can be satisfied.
1520
 * In the end we try to use one of the two constraints to eliminate
1521
 * a column.
1522
 *
1523
 * This function assumes that at least two more rows and at least
1524
 * two more elements in the constraint array are available in the tableau.
1525
 */
1526
static int add_lexmin_eq(struct isl_tab *tab, isl_int *eq) WARN_UNUSED;
1527
static int add_lexmin_eq(struct isl_tab *tab, isl_int *eq)
1528
301
{
1529
301
  int r1, r2;
1530
301
  int row;
1531
301
  struct isl_tab_undo *snap;
1532
301
1533
301
  if (!tab)
1534
0
    return -1;
1535
301
  snap = isl_tab_snap(tab);
1536
301
  r1 = isl_tab_add_row(tab, eq);
1537
301
  if (r1 < 0)
1538
0
    return -1;
1539
301
  tab->con[r1].is_nonneg = 1;
1540
301
  if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r1]) < 0)
1541
0
    return -1;
1542
301
1543
301
  row = tab->con[r1].index;
1544
301
  if (
is_constant(tab, row)301
) {
1545
54
    if (
!54
isl_int_is_zero54
(tab->mat->row[row][1]) ||
1546
54
        
(tab->M && 54
!0
isl_int_is_zero0
(tab->mat->row[row][2]))) {
1547
0
      if (isl_tab_mark_empty(tab) < 0)
1548
0
        return -1;
1549
0
      return 0;
1550
0
    }
1551
54
    
if (54
isl_tab_rollback(tab, snap) < 054
)
1552
0
      return -1;
1553
54
    return 0;
1554
54
  }
1555
247
1556
247
  
if (247
restore_lexmin(tab) < 0247
)
1557
0
    return -1;
1558
247
  
if (247
tab->empty247
)
1559
14
    return 0;
1560
233
1561
233
  isl_seq_neg(eq, eq, 1 + tab->n_var);
1562
233
1563
233
  r2 = isl_tab_add_row(tab, eq);
1564
233
  if (r2 < 0)
1565
0
    return -1;
1566
233
  tab->con[r2].is_nonneg = 1;
1567
233
  if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r2]) < 0)
1568
0
    return -1;
1569
233
1570
233
  
if (233
restore_lexmin(tab) < 0233
)
1571
0
    return -1;
1572
233
  
if (233
tab->empty233
)
1573
0
    return 0;
1574
233
1575
233
  
if (233
!tab->con[r1].is_row233
) {
1576
0
    if (isl_tab_kill_col(tab, tab->con[r1].index) < 0)
1577
0
      return -1;
1578
233
  } else 
if (233
!tab->con[r2].is_row233
) {
1579
0
    if (isl_tab_kill_col(tab, tab->con[r2].index) < 0)
1580
0
      return -1;
1581
233
  }
1582
233
1583
233
  
if (233
tab->bmap233
) {
1584
0
    tab->bmap = isl_basic_map_add_ineq(tab->bmap, eq);
1585
0
    if (isl_tab_push(tab, isl_tab_undo_bmap_ineq) < 0)
1586
0
      return -1;
1587
0
    isl_seq_neg(eq, eq, 1 + tab->n_var);
1588
0
    tab->bmap = isl_basic_map_add_ineq(tab->bmap, eq);
1589
0
    isl_seq_neg(eq, eq, 1 + tab->n_var);
1590
0
    if (isl_tab_push(tab, isl_tab_undo_bmap_ineq) < 0)
1591
0
      return -1;
1592
0
    
if (0
!tab->bmap0
)
1593
0
      return -1;
1594
233
  }
1595
233
1596
233
  return 0;
1597
233
}
1598
1599
/* Add an inequality to the tableau, resolving violations using
1600
 * restore_lexmin.
1601
 *
1602
 * This function assumes that at least one more row and at least
1603
 * one more element in the constraint array are available in the tableau.
1604
 */
1605
static struct isl_tab *add_lexmin_ineq(struct isl_tab *tab, isl_int *ineq)
1606
21.3k
{
1607
21.3k
  int r;
1608
21.3k
1609
21.3k
  if (!tab)
1610
0
    return NULL;
1611
21.3k
  
if (21.3k
tab->bmap21.3k
) {
1612
0
    tab->bmap = isl_basic_map_add_ineq(tab->bmap, ineq);
1613
0
    if (isl_tab_push(tab, isl_tab_undo_bmap_ineq) < 0)
1614
0
      goto error;
1615
0
    
if (0
!tab->bmap0
)
1616
0
      goto error;
1617
21.3k
  }
1618
21.3k
  r = isl_tab_add_row(tab, ineq);
1619
21.3k
  if (r < 0)
1620
0
    goto error;
1621
21.3k
  tab->con[r].is_nonneg = 1;
1622
21.3k
  if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r]) < 0)
1623
0
    goto error;
1624
21.3k
  
if (21.3k
isl_tab_row_is_redundant(tab, tab->con[r].index)21.3k
) {
1625
41
    if (isl_tab_mark_redundant(tab, tab->con[r].index) < 0)
1626
0
      goto error;
1627
41
    return tab;
1628
41
  }
1629
21.3k
1630
21.3k
  
if (21.3k
restore_lexmin(tab) < 021.3k
)
1631
0
    goto error;
1632
21.3k
  
if (21.3k
!tab->empty && 21.3k
tab->con[r].is_row20.8k
&&
1633
16.2k
     isl_tab_row_is_redundant(tab, tab->con[r].index))
1634
0
    
if (0
isl_tab_mark_redundant(tab, tab->con[r].index) < 00
)
1635
0
      goto error;
1636
21.3k
  return tab;
1637
0
error:
1638
0
  isl_tab_free(tab);
1639
0
  return NULL;
1640
21.3k
}
1641
1642
/* Check if the coefficients of the parameters are all integral.
1643
 */
1644
static int integer_parameter(struct isl_tab *tab, int row)
1645
21.1k
{
1646
21.1k
  int i;
1647
21.1k
  int col;
1648
21.1k
  unsigned off = 2 + tab->M;
1649
21.1k
1650
93.2k
  for (i = 0; 
i < tab->n_param93.2k
;
++i72.0k
) {
1651
72.5k
    /* Eliminated parameter */
1652
72.5k
    if (tab->var[i].is_row)
1653
36.5k
      continue;
1654
36.0k
    col = tab->var[i].index;
1655
36.0k
    if (
!36.0k
isl_int_is_divisible_by36.0k
(tab->mat->row[row][off + col],
1656
36.0k
            tab->mat->row[row][0]))
1657
519
      return 0;
1658
72.5k
  }
1659
26.0k
  
for (i = 0; 20.6k
i < tab->n_div26.0k
;
++i5.41k
) {
1660
5.48k
    if (tab->var[tab->n_var - tab->n_div + i].is_row)
1661
1.37k
      continue;
1662
4.10k
    col = tab->var[tab->n_var - tab->n_div + i].index;
1663
4.10k
    if (
!4.10k
isl_int_is_divisible_by4.10k
(tab->mat->row[row][off + col],
1664
4.10k
            tab->mat->row[row][0]))
1665
73
      return 0;
1666
5.48k
  }
1667
20.6k
  return 1;
1668
21.1k
}
1669
1670
/* Check if the coefficients of the non-parameter variables are all integral.
1671
 */
1672
static int integer_variable(struct isl_tab *tab, int row)
1673
1.07k
{
1674
1.07k
  int i;
1675
1.07k
  unsigned off = 2 + tab->M;
1676
1.07k
1677
2.99k
  for (i = tab->n_dead; 
i < tab->n_col2.99k
;
++i1.92k
) {
1678
2.90k
    if (tab->col_var[i] >= 0 &&
1679
1.73k
        (tab->col_var[i] < tab->n_param ||
1680
142
         tab->col_var[i] >= tab->n_var - tab->n_div))
1681
1.70k
      continue;
1682
1.20k
    
if (1.20k
!1.20k
isl_int_is_divisible_by1.20k
(tab->mat->row[row][off + i],
1683
1.20k
            tab->mat->row[row][0]))
1684
986
      return 0;
1685
2.90k
  }
1686
88
  return 1;
1687
1.07k
}
1688
1689
/* Check if the constant term is integral.
1690
 */
1691
static int integer_constant(struct isl_tab *tab, int row)
1692
21.1k
{
1693
21.1k
  return isl_int_is_divisible_by(tab->mat->row[row][1],
1694
21.1k
          tab->mat->row[row][0]);
1695
21.1k
}
1696
1697
#define I_CST 1 << 0
1698
#define I_PAR 1 << 1
1699
#define I_VAR 1 << 2
1700
1701
/* Check for next (non-parameter) variable after "var" (first if var == -1)
1702
 * that is non-integer and therefore requires a cut and return
1703
 * the index of the variable.
1704
 * For parametric tableaus, there are three parts in a row,
1705
 * the constant, the coefficients of the parameters and the rest.
1706
 * For each part, we check whether the coefficients in that part
1707
 * are all integral and if so, set the corresponding flag in *f.
1708
 * If the constant and the parameter part are integral, then the
1709
 * current sample value is integral and no cut is required
1710
 * (irrespective of whether the variable part is integral).
1711
 */
1712
static int next_non_integer_var(struct isl_tab *tab, int var, int *f)
1713
8.55k
{
1714
8.55k
  var = var < 0 ? 
tab->n_param8.55k
:
var + 10
;
1715
8.55k
1716
35.0k
  for (; 
var < tab->n_var - tab->n_div35.0k
;
++var26.5k
) {
1717
27.5k
    int flags = 0;
1718
27.5k
    int row;
1719
27.5k
    if (!tab->var[var].is_row)
1720
6.38k
      continue;
1721
21.1k
    row = tab->var[var].index;
1722
21.1k
    if (integer_constant(tab, row))
1723
20.4k
      ISL_FL_SET(flags, I_CST);
1724
21.1k
    if (integer_parameter(tab, row))
1725
20.6k
      ISL_FL_SET(flags, I_PAR);
1726
21.1k
    if (
ISL_FL_ISSET21.1k
(flags, I_CST) && 21.1k
ISL_FL_ISSET20.4k
(flags, I_PAR))
1727
20.1k
      continue;
1728
1.07k
    
if (1.07k
integer_variable(tab, row)1.07k
)
1729
88
      ISL_FL_SET(flags, I_VAR);
1730
27.5k
    *f = flags;
1731
27.5k
    return var;
1732
27.5k
  }
1733
7.47k
  return -1;
1734
8.55k
}
1735
1736
/* Check for first (non-parameter) variable that is non-integer and
1737
 * therefore requires a cut and return the corresponding row.
1738
 * For parametric tableaus, there are three parts in a row,
1739
 * the constant, the coefficients of the parameters and the rest.
1740
 * For each part, we check whether the coefficients in that part
1741
 * are all integral and if so, set the corresponding flag in *f.
1742
 * If the constant and the parameter part are integral, then the
1743
 * current sample value is integral and no cut is required
1744
 * (irrespective of whether the variable part is integral).
1745
 */
1746
static int first_non_integer_row(struct isl_tab *tab, int *f)
1747
7.77k
{
1748
7.77k
  int var = next_non_integer_var(tab, -1, f);
1749
7.77k
1750
7.77k
  return var < 0 ? 
-16.75k
:
tab->var[var].index1.01k
;
1751
7.77k
}
1752
1753
/* Add a (non-parametric) cut to cut away the non-integral sample
1754
 * value of the given row.
1755
 *
1756
 * If the row is given by
1757
 *
1758
 *  m r = f + \sum_i a_i y_i
1759
 *
1760
 * then the cut is
1761
 *
1762
 *  c = - {-f/m} + \sum_i {a_i/m} y_i >= 0
1763
 *
1764
 * The big parameter, if any, is ignored, since it is assumed to be big
1765
 * enough to be divisible by any integer.
1766
 * If the tableau is actually a parametric tableau, then this function
1767
 * is only called when all coefficients of the parameters are integral.
1768
 * The cut therefore has zero coefficients for the parameters.
1769
 *
1770
 * The current value is known to be negative, so row_sign, if it
1771
 * exists, is set accordingly.
1772
 *
1773
 * Return the row of the cut or -1.
1774
 */
1775
static int add_cut(struct isl_tab *tab, int row)
1776
482
{
1777
482
  int i;
1778
482
  int r;
1779
482
  isl_int *r_row;
1780
482
  unsigned off = 2 + tab->M;
1781
482
1782
482
  if (isl_tab_extend_cons(tab, 1) < 0)
1783
0
    return -1;
1784
482
  r = isl_tab_allocate_con(tab);
1785
482
  if (r < 0)
1786
0
    return -1;
1787
482
1788
482
  r_row = tab->mat->row[tab->con[r].index];
1789
482
  isl_int_set(r_row[0], tab->mat->row[row][0]);
1790
482
  isl_int_neg(r_row[1], tab->mat->row[row][1]);
1791
482
  isl_int_fdiv_r(r_row[1], r_row[1], tab->mat->row[row][0]);
1792
482
  isl_int_neg(r_row[1], r_row[1]);
1793
482
  if (tab->M)
1794
427
    isl_int_set_si(r_row[2], 0);
1795
4.51k
  for (i = 0; 
i < tab->n_col4.51k
;
++i4.03k
)
1796
4.03k
    isl_int_fdiv_r(r_row[off + i],
1797
482
      tab->mat->row[row][off + i], tab->mat->row[row][0]);
1798
482
1799
482
  tab->con[r].is_nonneg = 1;
1800
482
  if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r]) < 0)
1801
0
    return -1;
1802
482
  
if (482
tab->row_sign482
)
1803
427
    tab->row_sign[tab->con[r].index] = isl_tab_row_neg;
1804
482
1805
482
  return tab->con[r].index;
1806
482
}
1807
1808
0
#define CUT_ALL 1
1809
1.24k
#define CUT_ONE 0
1810
1811
/* Given a non-parametric tableau, add cuts until an integer
1812
 * sample point is obtained or until the tableau is determined
1813
 * to be integer infeasible.
1814
 * As long as there is any non-integer value in the sample point,
1815
 * we add appropriate cuts, if possible, for each of these
1816
 * non-integer values and then resolve the violated
1817
 * cut constraints using restore_lexmin.
1818
 * If one of the corresponding rows is equal to an integral
1819
 * combination of variables/constraints plus a non-integral constant,
1820
 * then there is no way to obtain an integer point and we return
1821
 * a tableau that is marked empty.
1822
 * The parameter cutting_strategy controls the strategy used when adding cuts
1823
 * to remove non-integer points. CUT_ALL adds all possible cuts
1824
 * before continuing the search. CUT_ONE adds only one cut at a time.
1825
 */
1826
static struct isl_tab *cut_to_integer_lexmin(struct isl_tab *tab,
1827
  int cutting_strategy)
1828
1.18k
{
1829
1.18k
  int var;
1830
1.18k
  int row;
1831
1.18k
  int flags;
1832
1.18k
1833
1.18k
  if (!tab)
1834
0
    return NULL;
1835
1.18k
  
if (1.18k
tab->empty1.18k
)
1836
463
    return tab;
1837
725
1838
778
  
while (725
(var = next_non_integer_var(tab, -1, &flags)) != -1778
) {
1839
55
    do {
1840
55
      if (
ISL_FL_ISSET55
(flags, I_VAR)) {
1841
0
        if (isl_tab_mark_empty(tab) < 0)
1842
0
          goto error;
1843
0
        return tab;
1844
0
      }
1845
55
      row = tab->var[var].index;
1846
55
      row = add_cut(tab, row);
1847
55
      if (row < 0)
1848
0
        goto error;
1849
55
      
if (55
cutting_strategy == 55
CUT_ONE55
)
1850
55
        break;
1851
0
    } while ((var = next_non_integer_var(tab, var, &flags)) != -1);
1852
55
    
if (55
restore_lexmin(tab) < 055
)
1853
0
      goto error;
1854
55
    
if (55
tab->empty55
)
1855
2
      break;
1856
55
  }
1857
725
  return tab;
1858
0
error:
1859
0
  isl_tab_free(tab);
1860
0
  return NULL;
1861
1.18k
}
1862
1863
/* Check whether all the currently active samples also satisfy the inequality
1864
 * "ineq" (treated as an equality if eq is set).
1865
 * Remove those samples that do not.
1866
 */
1867
static struct isl_tab *check_samples(struct isl_tab *tab, isl_int *ineq, int eq)
1868
14.8k
{
1869
14.8k
  int i;
1870
14.8k
  isl_int v;
1871
14.8k
1872
14.8k
  if (!tab)
1873
0
    return NULL;
1874
14.8k
1875
14.8k
  
isl_assert14.8k
(tab->mat->ctx, tab->bmap, goto error);
1876
14.8k
  
isl_assert14.8k
(tab->mat->ctx, tab->samples, goto error);
1877
14.8k
  
isl_assert14.8k
(tab->mat->ctx, tab->samples->n_col == 1 + tab->n_var, goto error);
1878
14.8k
1879
14.8k
  
isl_int_init14.8k
(v);
1880
37.0k
  for (i = tab->n_outside; 
i < tab->n_sample37.0k
;
++i22.1k
) {
1881
22.1k
    int sgn;
1882
22.1k
    isl_seq_inner_product(ineq, tab->samples->row[i],
1883
22.1k
          1 + tab->n_var, &v);
1884
22.1k
    sgn = isl_int_sgn(v);
1885
22.1k
    if (
eq ? 22.1k
(sgn == 0)9.23k
:
(sgn >= 0)12.8k
)
1886
15.6k
      continue;
1887
6.50k
    tab = isl_tab_drop_sample(tab, i);
1888
6.50k
    if (!tab)
1889
0
      break;
1890
22.1k
  }
1891
14.8k
  isl_int_clear(v);
1892
14.8k
1893
14.8k
  return tab;
1894
0
error:
1895
0
  isl_tab_free(tab);
1896
0
  return NULL;
1897
14.8k
}
1898
1899
/* Check whether the sample value of the tableau is finite,
1900
 * i.e., either the tableau does not use a big parameter, or
1901
 * all values of the variables are equal to the big parameter plus
1902
 * some constant.  This constant is the actual sample value.
1903
 */
1904
static int sample_is_finite(struct isl_tab *tab)
1905
0
{
1906
0
  int i;
1907
0
1908
0
  if (!tab->M)
1909
0
    return 1;
1910
0
1911
0
  
for (i = 0; 0
i < tab->n_var0
;
++i0
) {
1912
0
    int row;
1913
0
    if (!tab->var[i].is_row)
1914
0
      return 0;
1915
0
    row = tab->var[i].index;
1916
0
    if (isl_int_ne(tab->mat->row[row][0], tab->mat->row[row][2]))
1917
0
      return 0;
1918
0
  }
1919
0
  return 1;
1920
0
}
1921
1922
/* Check if the context tableau of sol has any integer points.
1923
 * Leave tab in empty state if no integer point can be found.
1924
 * If an integer point can be found and if moreover it is finite,
1925
 * then it is added to the list of sample values.
1926
 *
1927
 * This function is only called when none of the currently active sample
1928
 * values satisfies the most recently added constraint.
1929
 */
1930
static struct isl_tab *check_integer_feasible(struct isl_tab *tab)
1931
0
{
1932
0
  struct isl_tab_undo *snap;
1933
0
1934
0
  if (!tab)
1935
0
    return NULL;
1936
0
1937
0
  snap = isl_tab_snap(tab);
1938
0
  if (isl_tab_push_basis(tab) < 0)
1939
0
    goto error;
1940
0
1941
0
  
tab = cut_to_integer_lexmin(tab, 0
CUT_ALL0
);
1942
0
  if (!tab)
1943
0
    goto error;
1944
0
1945
0
  
if (0
!tab->empty && 0
sample_is_finite(tab)0
) {
1946
0
    struct isl_vec *sample;
1947
0
1948
0
    sample = isl_tab_get_sample_value(tab);
1949
0
1950
0
    if (isl_tab_add_sample(tab, sample) < 0)
1951
0
      goto error;
1952
0
  }
1953
0
1954
0
  
if (0
!tab->empty && 0
isl_tab_rollback(tab, snap) < 00
)
1955
0
    goto error;
1956
0
1957
0
  return tab;
1958
0
error:
1959
0
  isl_tab_free(tab);
1960
0
  return NULL;
1961
0
}
1962
1963
/* Check if any of the currently active sample values satisfies
1964
 * the inequality "ineq" (an equality if eq is set).
1965
 */
1966
static int tab_has_valid_sample(struct isl_tab *tab, isl_int *ineq, int eq)
1967
27.5k
{
1968
27.5k
  int i;
1969
27.5k
  isl_int v;
1970
27.5k
1971
27.5k
  if (!tab)
1972
0
    return -1;
1973
27.5k
1974
27.5k
  
isl_assert27.5k
(tab->mat->ctx, tab->bmap, return -1);
1975
27.5k
  
isl_assert27.5k
(tab->mat->ctx, tab->samples, return -1);
1976
27.5k
  
isl_assert27.5k
(tab->mat->ctx, tab->samples->n_col == 1 + tab->n_var, return -1);
1977
27.5k
1978
27.5k
  
isl_int_init27.5k
(v);
1979
45.8k
  for (i = tab->n_outside; 
i < tab->n_sample45.8k
;
++i18.3k
) {
1980
27.5k
    int sgn;
1981
27.5k
    isl_seq_inner_product(ineq, tab->samples->row[i],
1982
27.5k
          1 + tab->n_var, &v);
1983
27.5k
    sgn = isl_int_sgn(v);
1984
27.5k
    if (
eq ? 27.5k
(sgn == 0)9.12k
:
(sgn >= 0)18.4k
)
1985
9.21k
      break;
1986
27.5k
  }
1987
27.5k
  isl_int_clear(v);
1988
27.5k
1989
27.5k
  return i < tab->n_sample;
1990
27.5k
}
1991
1992
/* Insert a div specified by "div" to the tableau "tab" at position "pos" and
1993
 * return isl_bool_true if the div is obviously non-negative.
1994
 */
1995
static isl_bool context_tab_insert_div(struct isl_tab *tab, int pos,
1996
  __isl_keep isl_vec *div,
1997
  isl_stat (*add_ineq)(void *user, isl_int *), void *user)
1998
565
{
1999
565
  int i;
2000
565
  int r;
2001
565
  struct isl_mat *samples;
2002
565
  int nonneg;
2003
565
2004
565
  r = isl_tab_insert_div(tab, pos, div, add_ineq, user);
2005
565
  if (r < 0)
2006
0
    return isl_bool_error;
2007
565
  nonneg = tab->var[r].is_nonneg;
2008
565
  tab->var[r].frozen = 1;
2009
565
2010
565
  samples = isl_mat_extend(tab->samples,
2011
565
      tab->n_sample, 1 + tab->n_var);
2012
565
  tab->samples = samples;
2013
565
  if (!samples)
2014
0
    return isl_bool_error;
2015
1.51k
  
for (i = tab->n_outside; 565
i < samples->n_row1.51k
;
++i953
) {
2016
953
    isl_seq_inner_product(div->el + 1, samples->row[i],
2017
953
      div->size - 1, &samples->row[i][samples->n_col - 1]);
2018
953
    isl_int_fdiv_q(samples->row[i][samples->n_col - 1],
2019
953
             samples->row[i][samples->n_col - 1], div->el[0]);
2020
953
  }
2021
565
  tab->samples = isl_mat_move_cols(tab->samples, 1 + pos,
2022
565
          1 + tab->n_var - 1, 1);
2023
565
  if (!tab->samples)
2024
0
    return isl_bool_error;
2025
565
2026
565
  return nonneg;
2027
565
}
2028
2029
/* Add a div specified by "div" to both the main tableau and
2030
 * the context tableau.  In case of the main tableau, we only
2031
 * need to add an extra div.  In the context tableau, we also
2032
 * need to express the meaning of the div.
2033
 * Return the index of the div or -1 if anything went wrong.
2034
 *
2035
 * The new integer division is added before any unknown integer
2036
 * divisions in the context to ensure that it does not get
2037
 * equated to some linear combination involving unknown integer
2038
 * divisions.
2039
 */
2040
static int add_div(struct isl_tab *tab, struct isl_context *context,
2041
  __isl_keep isl_vec *div)
2042
565
{
2043
565
  int r;
2044
565
  int pos;
2045
565
  isl_bool nonneg;
2046
565
  struct isl_tab *context_tab = context->op->peek_tab(context);
2047
565
2048
565
  if (
!tab || 565
!context_tab565
)
2049
0
    goto error;
2050
565
2051
565
  pos = context_tab->n_var - context->n_unknown;
2052
565
  if ((nonneg = context->op->insert_div(context, pos, div)) < 0)
2053
0
    goto error;
2054
565
2055
565
  
if (565
!context->op->is_ok(context)565
)
2056
0
    goto error;
2057
565
2058
565
  pos = tab->n_var - context->n_unknown;
2059
565
  if (isl_tab_extend_vars(tab, 1) < 0)
2060
0
    goto error;
2061
565
  r = isl_tab_insert_var(tab, pos);
2062
565
  if (r < 0)
2063
0
    goto error;
2064
565
  
if (565
nonneg565
)
2065
76
    tab->var[r].is_nonneg = 1;
2066
565
  tab->var[r].frozen = 1;
2067
565
  tab->n_div++;
2068
565
2069
565
  return tab->n_div - 1 - context->n_unknown;
2070
0
error:
2071
0
  context->op->invalidate(context);
2072
0
  return -1;
2073
565
}
2074
2075
static int find_div(struct isl_tab *tab, isl_int *div, isl_int denom)
2076
592
{
2077
592
  int i;
2078
592
  unsigned total = isl_basic_map_total_dim(tab->bmap);
2079
592
2080
1.28k
  for (i = 0; 
i < tab->bmap->n_div1.28k
;
++i697
) {
2081
724
    if (isl_int_ne(tab->bmap->div[i][0], denom))
2082
641
      continue;
2083
83
    
if (83
!isl_seq_eq(tab->bmap->div[i] + 1, div, 1 + total)83
)
2084
56
      continue;
2085
27
    return i;
2086
27
  }
2087
565
  return -1;
2088
592
}
2089
2090
/* Return the index of a div that corresponds to "div".
2091
 * We first check if we already have such a div and if not, we create one.
2092
 */
2093
static int get_div(struct isl_tab *tab, struct isl_context *context,
2094
  struct isl_vec *div)
2095
592
{
2096
592
  int d;
2097
592
  struct isl_tab *context_tab = context->op->peek_tab(context);
2098
592
2099
592
  if (!context_tab)
2100
0
    return -1;
2101
592
2102
592
  d = find_div(context_tab, div->el + 1, div->el[0]);
2103
592
  if (d != -1)
2104
27
    return d;
2105
565
2106
565
  return add_div(tab, context, div);
2107
565
}
2108
2109
/* Add a parametric cut to cut away the non-integral sample value
2110
 * of the give row.
2111
 * Let a_i be the coefficients of the constant term and the parameters
2112
 * and let b_i be the coefficients of the variables or constraints
2113
 * in basis of the tableau.
2114
 * Let q be the div q = floor(\sum_i {-a_i} y_i).
2115
 *
2116
 * The cut is expressed as
2117
 *
2118
 *  c = \sum_i -{-a_i} y_i + \sum_i {b_i} x_i + q >= 0
2119
 *
2120
 * If q did not already exist in the context tableau, then it is added first.
2121
 * If q is in a column of the main tableau then the "+ q" can be accomplished
2122
 * by setting the corresponding entry to the denominator of the constraint.
2123
 * If q happens to be in a row of the main tableau, then the corresponding
2124
 * row needs to be added instead (taking care of the denominators).
2125
 * Note that this is very unlikely, but perhaps not entirely impossible.
2126
 *
2127
 * The current value of the cut is known to be negative (or at least
2128
 * non-positive), so row_sign is set accordingly.
2129
 *
2130
 * Return the row of the cut or -1.
2131
 */
2132
static int add_parametric_cut(struct isl_tab *tab, int row,
2133
  struct isl_context *context)
2134
504
{
2135
504
  struct isl_vec *div;
2136
504
  int d;
2137
504
  int i;
2138
504
  int r;
2139
504
  isl_int *r_row;
2140
504
  int col;
2141
504
  int n;
2142
504
  unsigned off = 2 + tab->M;
2143
504
2144
504
  if (!context)
2145
0
    return -1;
2146
504
2147
504
  div = get_row_parameter_div(tab, row);
2148
504
  if (!div)
2149
0
    return -1;
2150
504
2151
504
  n = tab->n_div - context->n_unknown;
2152
504
  d = context->op->get_div(context, tab, div);
2153
504
  isl_vec_free(div);
2154
504
  if (d < 0)
2155
0
    return -1;
2156
504
2157
504
  
if (504
isl_tab_extend_cons(tab, 1) < 0504
)
2158
0
    return -1;
2159
504
  r = isl_tab_allocate_con(tab);
2160
504
  if (r < 0)
2161
0
    return -1;
2162
504
2163
504
  r_row = tab->mat->row[tab->con[r].index];
2164
504
  isl_int_set(r_row[0], tab->mat->row[row][0]);
2165
504
  isl_int_neg(r_row[1], tab->mat->row[row][1]);
2166
504
  isl_int_fdiv_r(r_row[1], r_row[1], tab->mat->row[row][0]);
2167
504
  isl_int_neg(r_row[1], r_row[1]);
2168
504
  if (tab->M)
2169
504
    isl_int_set_si(r_row[2], 0);
2170
1.79k
  for (i = 0; 
i < tab->n_param1.79k
;
++i1.28k
) {
2171
1.28k
    if (tab->var[i].is_row)
2172
63
      continue;
2173
1.22k
    col = tab->var[i].index;
2174
1.22k
    isl_int_neg(r_row[off + col], tab->mat->row[row][off + col]);
2175
1.22k
    isl_int_fdiv_r(r_row[off + col], r_row[off + col],
2176
1.22k
        tab->mat->row[row][0]);
2177
1.22k
    isl_int_neg(r_row[off + col], r_row[off + col]);
2178
1.28k
  }
2179
1.64k
  for (i = 0; 
i < tab->n_div1.64k
;
++i1.14k
) {
2180
1.14k
    if (tab->var[tab->n_var - tab->n_div + i].is_row)
2181
31
      continue;
2182
1.11k
    col = tab->var[tab->n_var - tab->n_div + i].index;
2183
1.11k
    isl_int_neg(r_row[off + col], tab->mat->row[row][off + col]);
2184
1.11k
    isl_int_fdiv_r(r_row[off + col], r_row[off + col],
2185
1.11k
        tab->mat->row[row][0]);
2186
1.11k
    isl_int_neg(r_row[off + col], r_row[off + col]);
2187
1.14k
  }
2188
3.90k
  for (i = 0; 
i < tab->n_col3.90k
;
++i3.40k
) {
2189
3.40k
    if (tab->col_var[i] >= 0 &&
2190
2.33k
        (tab->col_var[i] < tab->n_param ||
2191
1.11k
         tab->col_var[i] >= tab->n_var - tab->n_div))
2192
2.33k
      continue;
2193
1.07k
    
isl_int_fdiv_r1.07k
(r_row[off + i],
2194
1.07k
      tab->mat->row[row][off + i], tab->mat->row[row][0]);
2195
1.07k
  }
2196
504
  if (
tab->var[tab->n_var - tab->n_div + d].is_row504
) {
2197
1
    isl_int gcd;
2198
1
    int d_row = tab->var[tab->n_var - tab->n_div + d].index;
2199
1
    isl_int_init(gcd);
2200
1
    isl_int_gcd(gcd, tab->mat->row[d_row][0], r_row[0]);
2201
1
    isl_int_divexact(r_row[0], r_row[0], gcd);
2202
1
    isl_int_divexact(gcd, tab->mat->row[d_row][0], gcd);
2203
1
    isl_seq_combine(r_row + 1, gcd, r_row + 1,
2204
1
        r_row[0], tab->mat->row[d_row] + 1,
2205
1
        off - 1 + tab->n_col);
2206
1
    isl_int_mul(r_row[0], r_row[0], tab->mat->row[d_row][0]);
2207
1
    isl_int_clear(gcd);
2208
504
  } else {
2209
503
    col = tab->var[tab->n_var - tab->n_div + d].index;
2210
503
    isl_int_set(r_row[off + col], tab->mat->row[row][0]);
2211
503
  }
2212
504
2213
504
  tab->con[r].is_nonneg = 1;
2214
504
  if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r]) < 0)
2215
0
    return -1;
2216
504
  
if (504
tab->row_sign504
)
2217
504
    tab->row_sign[tab->con[r].index] = isl_tab_row_neg;
2218
504
2219
504
  row = tab->con[r].index;
2220
504
2221
504
  if (
d >= n && 504
context->op->detect_equalities(context, tab) < 0502
)
2222
0
    return -1;
2223
504
2224
504
  return row;
2225
504
}
2226
2227
/* Construct a tableau for bmap that can be used for computing
2228
 * the lexicographic minimum (or maximum) of bmap.
2229
 * If not NULL, then dom is the domain where the minimum
2230
 * should be computed.  In this case, we set up a parametric
2231
 * tableau with row signs (initialized to "unknown").
2232
 * If M is set, then the tableau will use a big parameter.
2233
 * If max is set, then a maximum should be computed instead of a minimum.
2234
 * This means that for each variable x, the tableau will contain the variable
2235
 * x' = M - x, rather than x' = M + x.  This in turn means that the coefficient
2236
 * of the variables in all constraints are negated prior to adding them
2237
 * to the tableau.
2238
 */
2239
static __isl_give struct isl_tab *tab_for_lexmin(__isl_keep isl_basic_map *bmap,
2240
  __isl_keep isl_basic_set *dom, unsigned M, int max)
2241
6.96k
{
2242
6.96k
  int i;
2243
6.96k
  struct isl_tab *tab;
2244
6.96k
  unsigned n_var;
2245
6.96k
  unsigned o_var;
2246
6.96k
2247
6.96k
  tab = isl_tab_alloc(bmap->ctx, 2 * bmap->n_eq + bmap->n_ineq + 1,
2248
6.96k
          isl_basic_map_total_dim(bmap), M);
2249
6.96k
  if (!tab)
2250
0
    return NULL;
2251
6.96k
2252
6.96k
  
tab->rational = 6.96k
ISL_F_ISSET6.96k
(bmap, ISL_BASIC_MAP_RATIONAL);
2253
6.96k
  if (
dom6.96k
) {
2254
6.52k
    tab->n_param = isl_basic_set_total_dim(dom) - dom->n_div;
2255
6.52k
    tab->n_div = dom->n_div;
2256
6.52k
    tab->row_sign = isl_calloc_array(bmap->ctx,
2257
6.52k
          enum isl_tab_row_sign, tab->mat->n_row);
2258
6.52k
    if (
tab->mat->n_row && 6.52k
!tab->row_sign6.52k
)
2259
0
      goto error;
2260
6.96k
  }
2261
6.96k
  
if (6.96k
ISL_F_ISSET6.96k
(bmap, ISL_BASIC_MAP_EMPTY)) {
2262
0
    if (isl_tab_mark_empty(tab) < 0)
2263
0
      goto error;
2264
0
    return tab;
2265
0
  }
2266
6.96k
2267
27.8k
  
for (i = tab->n_param; 6.96k
i < tab->n_var - tab->n_div27.8k
;
++i20.8k
) {
2268
20.8k
    tab->var[i].is_nonneg = 1;
2269
20.8k
    tab->var[i].frozen = 1;
2270
20.8k
  }
2271
6.96k
  o_var = 1 + tab->n_param;
2272
6.96k
  n_var = tab->n_var - tab->n_param - tab->n_div;
2273
28.7k
  for (i = 0; 
i < bmap->n_eq28.7k
;
++i21.7k
) {
2274
21.7k
    if (max)
2275
17.3k
      isl_seq_neg(bmap->eq[i] + o_var,
2276
17.3k
            bmap->eq[i] + o_var, n_var);
2277
21.7k
    tab = add_lexmin_valid_eq(tab, bmap->eq[i]);
2278
21.7k
    if (max)
2279
17.3k
      isl_seq_neg(bmap->eq[i] + o_var,
2280
17.3k
            bmap->eq[i] + o_var, n_var);
2281
21.7k
    if (
!tab || 21.7k
tab->empty21.7k
)
2282
0
      return tab;
2283
21.7k
  }
2284
6.96k
  
if (6.96k
bmap->n_eq && 6.96k
restore_lexmin(tab) < 05.97k
)
2285
0
    goto error;
2286
27.5k
  
for (i = 0; 6.96k
i < bmap->n_ineq27.5k
;
++i20.6k
) {
2287
20.6k
    if (max)
2288
13.0k
      isl_seq_neg(bmap->ineq[i] + o_var,
2289
13.0k
            bmap->ineq[i] + o_var, n_var);
2290
20.6k
    tab = add_lexmin_ineq(tab, bmap->ineq[i]);
2291
20.6k
    if (max)
2292
13.0k
      isl_seq_neg(bmap->ineq[i] + o_var,
2293
13.0k
            bmap->ineq[i] + o_var, n_var);
2294
20.6k
    if (
!tab || 20.6k
tab->empty20.6k
)
2295
0
      return tab;
2296
20.6k
  }
2297
6.96k
  return tab;
2298
0
error:
2299
0
  isl_tab_free(tab);
2300
0
  return NULL;
2301
6.96k
}
2302
2303
/* Given a main tableau where more than one row requires a split,
2304
 * determine and return the "best" row to split on.
2305
 *
2306
 * Given two rows in the main tableau, if the inequality corresponding
2307
 * to the first row is redundant with respect to that of the second row
2308
 * in the current tableau, then it is better to split on the second row,
2309
 * since in the positive part, both rows will be positive.
2310
 * (In the negative part a pivot will have to be performed and just about
2311
 * anything can happen to the sign of the other row.)
2312
 *
2313
 * As a simple heuristic, we therefore select the row that makes the most
2314
 * of the other rows redundant.
2315
 *
2316
 * Perhaps it would also be useful to look at the number of constraints
2317
 * that conflict with any given constraint.
2318
 *
2319
 * best is the best row so far (-1 when we have not found any row yet).
2320
 * best_r is the number of other rows made redundant by row best.
2321
 * When best is still -1, bset_r is meaningless, but it is initialized
2322
 * to some arbitrary value (0) anyway.  Without this redundant initialization
2323
 * valgrind may warn about uninitialized memory accesses when isl
2324
 * is compiled with some versions of gcc.
2325
 */
2326
static int best_split(struct isl_tab *tab, struct isl_tab *context_tab)
2327
331
{
2328
331
  struct isl_tab_undo *snap;
2329
331
  int split;
2330
331
  int row;
2331
331
  int best = -1;
2332
331
  int best_r = 0;
2333
331
2334
331
  if (isl_tab_extend_cons(context_tab, 2) < 0)
2335
0
    return -1;
2336
331
2337
331
  snap = isl_tab_snap(context_tab);
2338
331
2339
2.66k
  for (split = tab->n_redundant; 
split < tab->n_row2.66k
;
++split2.33k
) {
2340
2.33k
    struct isl_tab_undo *snap2;
2341
2.33k
    struct isl_vec *ineq = NULL;
2342
2.33k
    int r = 0;
2343
2.33k
    int ok;
2344
2.33k
2345
2.33k
    if (!isl_tab_var_from_row(tab, split)->is_nonneg)
2346
0
      continue;
2347
2.33k
    
if (2.33k
tab->row_sign[split] != isl_tab_row_any2.33k
)
2348
1.53k
      continue;
2349
792
2350
792
    ineq = get_row_parameter_ineq(tab, split);
2351
792
    if (!ineq)
2352
0
      return -1;
2353
792
    ok = isl_tab_add_ineq(context_tab, ineq->el) >= 0;
2354
792
    isl_vec_free(ineq);
2355
792
    if (!ok)
2356
0
      return -1;
2357
792
2358
792
    snap2 = isl_tab_snap(context_tab);
2359
792
2360
6.69k
    for (row = tab->n_redundant; 
row < tab->n_row6.69k
;
++row5.90k
) {
2361
5.90k
      struct isl_tab_var *var;
2362
5.90k
2363
5.90k
      if (row == split)
2364
792
        continue;
2365
5.11k
      
if (5.11k
!isl_tab_var_from_row(tab, row)->is_nonneg5.11k
)
2366
0
        continue;
2367
5.11k
      
if (5.11k
tab->row_sign[row] != isl_tab_row_any5.11k
)
2368
3.75k
        continue;
2369
1.36k
2370
1.36k
      ineq = get_row_parameter_ineq(tab, row);
2371
1.36k
      if (!ineq)
2372
0
        return -1;
2373
1.36k
      ok = isl_tab_add_ineq(context_tab, ineq->el) >= 0;
2374
1.36k
      isl_vec_free(ineq);
2375
1.36k
      if (!ok)
2376
0
        return -1;
2377
1.36k
      var = &context_tab->con[context_tab->n_con - 1];
2378
1.36k
      if (!context_tab->empty &&
2379
1.36k
          !isl_tab_min_at_most_neg_one(context_tab, var))
2380
23
        r++;
2381
1.36k
      if (isl_tab_rollback(context_tab, snap2) < 0)
2382
0
        return -1;
2383
5.90k
    }
2384
792
    
if (792
best == -1 || 792
r > best_r461
) {
2385
335
      best = split;
2386
335
      best_r = r;
2387
335
    }
2388
792
    if (isl_tab_rollback(context_tab, snap) < 0)
2389
0
      return -1;
2390
2.33k
  }
2391
331
2392
331
  return best;
2393
331
}
2394
2395
static struct isl_basic_set *context_lex_peek_basic_set(
2396
  struct isl_context *context)
2397
0
{
2398
0
  struct isl_context_lex *clex = (struct isl_context_lex *)context;
2399
0
  if (!clex->tab)
2400
0
    return NULL;
2401
0
  return isl_tab_peek_bset(clex->tab);
2402
0
}
2403
2404
static struct isl_tab *context_lex_peek_tab(struct isl_context *context)
2405
0
{
2406
0
  struct isl_context_lex *clex = (struct isl_context_lex *)context;
2407
0
  return clex->tab;
2408
0
}
2409
2410
static void context_lex_add_eq(struct isl_context *context, isl_int *eq,
2411
    int check, int update)
2412
0
{
2413
0
  struct isl_context_lex *clex = (struct isl_context_lex *)context;
2414
0
  if (isl_tab_extend_cons(clex->tab, 2) < 0)
2415
0
    goto error;
2416
0
  
if (0
add_lexmin_eq(clex->tab, eq) < 00
)
2417
0
    goto error;
2418
0
  
if (0
check0
) {
2419
0
    int v = tab_has_valid_sample(clex->tab, eq, 1);
2420
0
    if (v < 0)
2421
0
      goto error;
2422
0
    
if (0
!v0
)
2423
0
      clex->tab = check_integer_feasible(clex->tab);
2424
0
  }
2425
0
  
if (0
update0
)
2426
0
    clex->tab = check_samples(clex->tab, eq, 1);
2427
0
  return;
2428
0
error:
2429
0
  isl_tab_free(clex->tab);
2430
0
  clex->tab = NULL;
2431
0
}
2432
2433
static void context_lex_add_ineq(struct isl_context *context, isl_int *ineq,
2434
    int check, int update)
2435
0
{
2436
0
  struct isl_context_lex *clex = (struct isl_context_lex *)context;
2437
0
  if (isl_tab_extend_cons(clex->tab, 1) < 0)
2438
0
    goto error;
2439
0
  clex->tab = add_lexmin_ineq(clex->tab, ineq);
2440
0
  if (
check0
) {
2441
0
    int v = tab_has_valid_sample(clex->tab, ineq, 0);
2442
0
    if (v < 0)
2443
0
      goto error;
2444
0
    
if (0
!v0
)
2445
0
      clex->tab = check_integer_feasible(clex->tab);
2446
0
  }
2447
0
  
if (0
update0
)
2448
0
    clex->tab = check_samples(clex->tab, ineq, 0);
2449
0
  return;
2450
0
error:
2451
0
  isl_tab_free(clex->tab);
2452
0
  clex->tab = NULL;
2453
0
}
2454
2455
static isl_stat context_lex_add_ineq_wrap(void *user, isl_int *ineq)
2456
0
{
2457
0
  struct isl_context *context = (struct isl_context *)user;
2458
0
  context_lex_add_ineq(context, ineq, 0, 0);
2459
0
  return context->op->is_ok(context) ? 
isl_stat_ok0
:
isl_stat_error0
;
2460
0
}
2461
2462
/* Check which signs can be obtained by "ineq" on all the currently
2463
 * active sample values.  See row_sign for more information.
2464
 */
2465
static enum isl_tab_row_sign tab_ineq_sign(struct isl_tab *tab, isl_int *ineq,
2466
  int strict)
2467
11.3k
{
2468
11.3k
  int i;
2469
11.3k
  int sgn;
2470
11.3k
  isl_int tmp;
2471
11.3k
  enum isl_tab_row_sign res = isl_tab_row_unknown;
2472
11.3k
2473
11.3k
  isl_assert(tab->mat->ctx, tab->samples, return isl_tab_row_unknown);
2474
11.3k
  
isl_assert11.3k
(tab->mat->ctx, tab->samples->n_col == 1 + tab->n_var,
2475
11.3k
      return isl_tab_row_unknown);
2476
11.3k
2477
11.3k
  
isl_int_init11.3k
(tmp);
2478
26.9k
  for (i = tab->n_outside; 
i < tab->n_sample26.9k
;
++i15.6k
) {
2479
16.1k
    isl_seq_inner_product(tab->samples->row[i], ineq,
2480
16.1k
          1 + tab->n_var, &tmp);
2481
16.1k
    sgn = isl_int_sgn(tmp);
2482
16.1k
    if (
sgn > 0 || 16.1k
(sgn == 0 && 7.20k
strict4.26k
)) {
2483
12.8k
      if (res == isl_tab_row_unknown)
2484
8.87k
        res = isl_tab_row_pos;
2485
12.8k
      if (res == isl_tab_row_neg)
2486
244
        res = isl_tab_row_any;
2487
12.8k
    }
2488
16.1k
    if (
sgn < 016.1k
) {
2489
2.93k
      if (res == isl_tab_row_unknown)
2490
2.37k
        res = isl_tab_row_neg;
2491
2.93k
      if (res == isl_tab_row_pos)
2492
292
        res = isl_tab_row_any;
2493
2.93k
    }
2494
16.1k
    if (res == isl_tab_row_any)
2495
536
      break;
2496
16.1k
  }
2497
11.3k
  isl_int_clear(tmp);
2498
11.3k
2499
11.3k
  return res;
2500
11.3k
}
2501
2502
static enum isl_tab_row_sign context_lex_ineq_sign(struct isl_context *context,
2503
      isl_int *ineq, int strict)
2504
0
{
2505
0
  struct isl_context_lex *clex = (struct isl_context_lex *)context;
2506
0
  return tab_ineq_sign(clex->tab, ineq, strict);
2507
0
}
2508
2509
/* Check whether "ineq" can be added to the tableau without rendering
2510
 * it infeasible.
2511
 */
2512
static int context_lex_test_ineq(struct isl_context *context, isl_int *ineq)
2513
0
{
2514
0
  struct isl_context_lex *clex = (struct isl_context_lex *)context;
2515
0
  struct isl_tab_undo *snap;
2516
0
  int feasible;
2517
0
2518
0
  if (!clex->tab)
2519
0
    return -1;
2520
0
2521
0
  
if (0
isl_tab_extend_cons(clex->tab, 1) < 00
)
2522
0
    return -1;
2523
0
2524
0
  snap = isl_tab_snap(clex->tab);
2525
0
  if (isl_tab_push_basis(clex->tab) < 0)
2526
0
    return -1;
2527
0
  clex->tab = add_lexmin_ineq(clex->tab, ineq);
2528
0
  clex->tab = check_integer_feasible(clex->tab);
2529
0
  if (!clex->tab)
2530
0
    return -1;
2531
0
  feasible = !clex->tab->empty;
2532
0
  if (isl_tab_rollback(clex->tab, snap) < 0)
2533
0
    return -1;
2534
0
2535
0
  return feasible;
2536
0
}
2537
2538
static int context_lex_get_div(struct isl_context *context, struct isl_tab *tab,
2539
    struct isl_vec *div)
2540
0
{
2541
0
  return get_div(tab, context, div);
2542
0
}
2543
2544
/* Insert a div specified by "div" to the context tableau at position "pos" and
2545
 * return isl_bool_true if the div is obviously non-negative.
2546
 * context_tab_add_div will always return isl_bool_true, because all variables
2547
 * in a isl_context_lex tableau are non-negative.
2548
 * However, if we are using a big parameter in the context, then this only
2549
 * reflects the non-negativity of the variable used to _encode_ the
2550
 * div, i.e., div' = M + div, so we can't draw any conclusions.
2551
 */
2552
static isl_bool context_lex_insert_div(struct isl_context *context, int pos,
2553
  __isl_keep isl_vec *div)
2554
0
{
2555
0
  struct isl_context_lex *clex = (struct isl_context_lex *)context;
2556
0
  isl_bool nonneg;
2557
0
  nonneg = context_tab_insert_div(clex->tab, pos, div,
2558
0
          context_lex_add_ineq_wrap, context);
2559
0
  if (nonneg < 0)
2560
0
    return isl_bool_error;
2561
0
  
if (0
clex->tab->M0
)
2562
0
    return isl_bool_false;
2563
0
  return nonneg;
2564
0
}
2565
2566
static int context_lex_detect_equalities(struct isl_context *context,
2567
    struct isl_tab *tab)
2568
0
{
2569
0
  return 0;
2570
0
}
2571
2572
static int context_lex_best_split(struct isl_context *context,
2573
    struct isl_tab *tab)
2574
0
{
2575
0
  struct isl_context_lex *clex = (struct isl_context_lex *)context;
2576
0
  struct isl_tab_undo *snap;
2577
0
  int r;
2578
0
2579
0
  snap = isl_tab_snap(clex->tab);
2580
0
  if (isl_tab_push_basis(clex->tab) < 0)
2581
0
    return -1;
2582
0
  r = best_split(tab, clex->tab);
2583
0
2584
0
  if (
r >= 0 && 0
isl_tab_rollback(clex->tab, snap) < 00
)
2585
0
    return -1;
2586
0
2587
0
  return r;
2588
0
}
2589
2590
static int context_lex_is_empty(struct isl_context *context)
2591
0
{
2592
0
  struct isl_context_lex *clex = (struct isl_context_lex *)context;
2593
0
  if (!clex->tab)
2594
0
    return -1;
2595
0
  return clex->tab->empty;
2596
0
}
2597
2598
static void *context_lex_save(struct isl_context *context)
2599
0
{
2600
0
  struct isl_context_lex *clex = (struct isl_context_lex *)context;
2601
0
  struct isl_tab_undo *snap;
2602
0
2603
0
  snap = isl_tab_snap(clex->tab);
2604
0
  if (isl_tab_push_basis(clex->tab) < 0)
2605
0
    return NULL;
2606
0
  
if (0
isl_tab_save_samples(clex->tab) < 00
)
2607
0
    return NULL;
2608
0
2609
0
  return snap;
2610
0
}
2611
2612
static void context_lex_restore(struct isl_context *context, void *save)
2613
0
{
2614
0
  struct isl_context_lex *clex = (struct isl_context_lex *)context;
2615
0
  if (
isl_tab_rollback(clex->tab, (struct isl_tab_undo *)save) < 00
) {
2616
0
    isl_tab_free(clex->tab);
2617
0
    clex->tab = NULL;
2618
0
  }
2619
0
}
2620
2621
static void context_lex_discard(void *save)
2622
0
{
2623
0
}
2624
2625
static int context_lex_is_ok(struct isl_context *context)
2626
0
{
2627
0
  struct isl_context_lex *clex = (struct isl_context_lex *)context;
2628
0
  return !!clex->tab;
2629
0
}
2630
2631
/* For each variable in the context tableau, check if the variable can
2632
 * only attain non-negative values.  If so, mark the parameter as non-negative
2633
 * in the main tableau.  This allows for a more direct identification of some
2634
 * cases of violated constraints.
2635
 */
2636
static struct isl_tab *tab_detect_nonnegative_parameters(struct isl_tab *tab,
2637
  struct isl_tab *context_tab)
2638
6.52k
{
2639
6.52k
  int i;
2640
6.52k
  struct isl_tab_undo *snap;
2641
6.52k
  struct isl_vec *ineq = NULL;
2642
6.52k
  struct isl_tab_var *var;
2643
6.52k
  int n;
2644
6.52k
2645
6.52k
  if (context_tab->n_var == 0)
2646
1.40k
    return tab;
2647
5.12k
2648
5.12k
  ineq = isl_vec_alloc(tab->mat->ctx, 1 + context_tab->n_var);
2649
5.12k
  if (!ineq)
2650
0
    goto error;
2651
5.12k
2652
5.12k
  
if (5.12k
isl_tab_extend_cons(context_tab, 1) < 05.12k
)
2653
0
    goto error;
2654
5.12k
2655
5.12k
  snap = isl_tab_snap(context_tab);
2656
5.12k
2657
5.12k
  n = 0;
2658
5.12k
  isl_seq_clr(ineq->el, ineq->size);
2659
25.0k
  for (i = 0; 
i < context_tab->n_var25.0k
;
++i19.9k
) {
2660
19.9k
    isl_int_set_si(ineq->el[1 + i], 1);
2661
19.9k
    if (isl_tab_add_ineq(context_tab, ineq->el) < 0)
2662
0
      goto error;
2663
19.9k
    var = &context_tab->con[context_tab->n_con - 1];
2664
19.9k
    if (!context_tab->empty &&
2665
19.9k
        
!isl_tab_min_at_most_neg_one(context_tab, var)19.6k
) {
2666
14.7k
      int j = i;
2667
14.7k
      if (i >= tab->n_param)
2668
95
        j = i - tab->n_param + tab->n_var - tab->n_div;
2669
14.7k
      tab->var[j].is_nonneg = 1;
2670
14.7k
      n++;
2671
14.7k
    }
2672
19.9k
    isl_int_set_si(ineq->el[1 + i], 0);
2673
19.9k
    if (isl_tab_rollback(context_tab, snap) < 0)
2674
0
      goto error;
2675
19.9k
  }
2676
5.12k
2677
5.12k
  
if (5.12k
context_tab->M && 5.12k
n == context_tab->n_var0
) {
2678
0
    context_tab->mat = isl_mat_drop_cols(context_tab->mat, 2, 1);
2679
0
    context_tab->M = 0;
2680
0
  }
2681
5.12k
2682
5.12k
  isl_vec_free(ineq);
2683
5.12k
  return tab;
2684
0
error:
2685
0
  isl_vec_free(ineq);
2686
0
  isl_tab_free(tab);
2687
0
  return NULL;
2688
6.52k
}
2689
2690
static struct isl_tab *context_lex_detect_nonnegative_parameters(
2691
  struct isl_context *context, struct isl_tab *tab)
2692
0
{
2693
0
  struct isl_context_lex *clex = (struct isl_context_lex *)context;
2694
0
  struct isl_tab_undo *snap;
2695
0
2696
0
  if (!tab)
2697
0
    return NULL;
2698
0
2699
0
  snap = isl_tab_snap(clex->tab);
2700
0
  if (isl_tab_push_basis(clex->tab) < 0)
2701
0
    goto error;
2702
0
2703
0
  tab = tab_detect_nonnegative_parameters(tab, clex->tab);
2704
0
2705
0
  if (isl_tab_rollback(clex->tab, snap) < 0)
2706
0
    goto error;
2707
0
2708
0
  return tab;
2709
0
error:
2710
0
  isl_tab_free(tab);
2711
0
  return NULL;
2712
0
}
2713
2714
static void context_lex_invalidate(struct isl_context *context)
2715
0
{
2716
0
  struct isl_context_lex *clex = (struct isl_context_lex *)context;
2717
0
  isl_tab_free(clex->tab);
2718
0
  clex->tab = NULL;
2719
0
}
2720
2721
static __isl_null struct isl_context *context_lex_free(
2722
  struct isl_context *context)
2723
0
{
2724
0
  struct isl_context_lex *clex = (struct isl_context_lex *)context;
2725
0
  isl_tab_free(clex->tab);
2726
0
  free(clex);
2727
0
2728
0
  return NULL;
2729
0
}
2730
2731
struct isl_context_op isl_context_lex_op = {
2732
  context_lex_detect_nonnegative_parameters,
2733
  context_lex_peek_basic_set,
2734
  context_lex_peek_tab,
2735
  context_lex_add_eq,
2736
  context_lex_add_ineq,
2737
  context_lex_ineq_sign,
2738
  context_lex_test_ineq,
2739
  context_lex_get_div,
2740
  context_lex_insert_div,
2741
  context_lex_detect_equalities,
2742
  context_lex_best_split,
2743
  context_lex_is_empty,
2744
  context_lex_is_ok,
2745
  context_lex_save,
2746
  context_lex_restore,
2747
  context_lex_discard,
2748
  context_lex_invalidate,
2749
  context_lex_free,
2750
};
2751
2752
static struct isl_tab *context_tab_for_lexmin(__isl_take isl_basic_set *bset)
2753
0
{
2754
0
  struct isl_tab *tab;
2755
0
2756
0
  if (!bset)
2757
0
    return NULL;
2758
0
  tab = tab_for_lexmin(bset_to_bmap(bset), NULL, 1, 0);
2759
0
  if (isl_tab_track_bset(tab, bset) < 0)
2760
0
    goto error;
2761
0
  tab = isl_tab_init_samples(tab);
2762
0
  return tab;
2763
0
error:
2764
0
  isl_tab_free(tab);
2765
0
  return NULL;
2766
0
}
2767
2768
static struct isl_context *isl_context_lex_alloc(struct isl_basic_set *dom)
2769
0
{
2770
0
  struct isl_context_lex *clex;
2771
0
2772
0
  if (!dom)
2773
0
    return NULL;
2774
0
2775
0
  
clex = 0
isl_alloc_type0
(dom->ctx, struct isl_context_lex);
2776
0
  if (!clex)
2777
0
    return NULL;
2778
0
2779
0
  clex->context.op = &isl_context_lex_op;
2780
0
2781
0
  clex->tab = context_tab_for_lexmin(isl_basic_set_copy(dom));
2782
0
  if (restore_lexmin(clex->tab) < 0)
2783
0
    goto error;
2784
0
  clex->tab = check_integer_feasible(clex->tab);
2785
0
  if (!clex->tab)
2786
0
    goto error;
2787
0
2788
0
  return &clex->context;
2789
0
error:
2790
0
  clex->context.op->free(&clex->context);
2791
0
  return NULL;
2792
0
}
2793
2794
/* Representation of the context when using generalized basis reduction.
2795
 *
2796
 * "shifted" contains the offsets of the unit hypercubes that lie inside the
2797
 * context.  Any rational point in "shifted" can therefore be rounded
2798
 * up to an integer point in the context.
2799
 * If the context is constrained by any equality, then "shifted" is not used
2800
 * as it would be empty.
2801
 */
2802
struct isl_context_gbr {
2803
  struct isl_context context;
2804
  struct isl_tab *tab;
2805
  struct isl_tab *shifted;
2806
  struct isl_tab *cone;
2807
};
2808
2809
static struct isl_tab *context_gbr_detect_nonnegative_parameters(
2810
  struct isl_context *context, struct isl_tab *tab)
2811
6.52k
{
2812
6.52k
  struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
2813
6.52k
  if (!tab)
2814
0
    return NULL;
2815
6.52k
  return tab_detect_nonnegative_parameters(tab, cgbr->tab);
2816
6.52k
}
2817
2818
static struct isl_basic_set *context_gbr_peek_basic_set(
2819
  struct isl_context *context)
2820
23.7k
{
2821
23.7k
  struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
2822
23.7k
  if (!cgbr->tab)
2823
0
    return NULL;
2824
23.7k
  return isl_tab_peek_bset(cgbr->tab);
2825
23.7k
}
2826
2827
static struct isl_tab *context_gbr_peek_tab(struct isl_context *context)
2828
31.5k
{
2829
31.5k
  struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
2830
31.5k
  return cgbr->tab;
2831
31.5k
}
2832
2833
/* Initialize the "shifted" tableau of the context, which
2834
 * contains the constraints of the original tableau shifted
2835
 * by the sum of all negative coefficients.  This ensures
2836
 * that any rational point in the shifted tableau can
2837
 * be rounded up to yield an integer point in the original tableau.
2838
 */
2839
static void gbr_init_shifted(struct isl_context_gbr *cgbr)
2840
154
{
2841
154
  int i, j;
2842
154
  struct isl_vec *cst;
2843
154
  struct isl_basic_set *bset = isl_tab_peek_bset(cgbr->tab);
2844
154
  unsigned dim = isl_basic_set_total_dim(bset);
2845
154
2846
154
  cst = isl_vec_alloc(cgbr->tab->mat->ctx, bset->n_ineq);
2847
154
  if (!cst)
2848
0
    return;
2849
154
2850
1.47k
  
for (i = 0; 154
i < bset->n_ineq1.47k
;
++i1.32k
) {
2851
1.32k
    isl_int_set(cst->el[i], bset->ineq[i][0]);
2852
10.1k
    for (j = 0; 
j < dim10.1k
;
++j8.77k
) {
2853
8.77k
      if (
!8.77k
isl_int_is_neg8.77k
(bset->ineq[i][1 + j]))
2854
7.57k
        continue;
2855
1.20k
      
isl_int_add1.20k
(bset->ineq[i][0], bset->ineq[i][0],
2856
1.20k
            bset->ineq[i][1 + j]);
2857
1.20k
    }
2858
1.32k
  }
2859
154
2860
154
  cgbr->shifted = isl_tab_from_basic_set(bset, 0);
2861
154
2862
1.47k
  for (i = 0; 
i < bset->n_ineq1.47k
;
++i1.32k
)
2863
1.32k
    isl_int_set(bset->ineq[i][0], cst->el[i]);
2864
154
2865
154
  isl_vec_free(cst);
2866
154
}
2867
2868
/* Check if the shifted tableau is non-empty, and if so
2869
 * use the sample point to construct an integer point
2870
 * of the context tableau.
2871
 */
2872
static struct isl_vec *gbr_get_shifted_sample(struct isl_context_gbr *cgbr)
2873
191
{
2874
191
  struct isl_vec *sample;
2875
191
2876
191
  if (!cgbr->shifted)
2877
154
    gbr_init_shifted(cgbr);
2878
191
  if (!cgbr->shifted)
2879
0
    return NULL;
2880
191
  
if (191
cgbr->shifted->empty191
)
2881
149
    return isl_vec_alloc(cgbr->tab->mat->ctx, 0);
2882
42
2883
42
  sample = isl_tab_get_sample_value(cgbr->shifted);
2884
42
  sample = isl_vec_ceil(sample);
2885
42
2886
42
  return sample;
2887
42
}
2888
2889
static __isl_give isl_basic_set *drop_constant_terms(
2890
  __isl_take isl_basic_set *bset)
2891
795
{
2892
795
  int i;
2893
795
2894
795
  if (!bset)
2895
0
    return NULL;
2896
795
2897
2.06k
  
for (i = 0; 795
i < bset->n_eq2.06k
;
++i1.26k
)
2898
1.26k
    isl_int_set_si(bset->eq[i][0], 0);
2899
795
2900
6.09k
  for (i = 0; 
i < bset->n_ineq6.09k
;
++i5.29k
)
2901
5.29k
    isl_int_set_si(bset->ineq[i][0], 0);
2902
795
2903
795
  return bset;
2904
795
}
2905
2906
static int use_shifted(struct isl_context_gbr *cgbr)
2907
1.00k
{
2908
1.00k
  if (!cgbr->tab)
2909
0
    return 0;
2910
1.00k
  
return cgbr->tab->bmap->n_eq == 0 && 1.00k
cgbr->tab->bmap->n_div == 0375
;
2911
1.00k
}
2912
2913
static struct isl_vec *gbr_get_sample(struct isl_context_gbr *cgbr)
2914
10.9k
{
2915
10.9k
  struct isl_basic_set *bset;
2916
10.9k
  struct isl_basic_set *cone;
2917
10.9k
2918
10.9k
  if (isl_tab_sample_is_integer(cgbr->tab))
2919
9.91k
    return isl_tab_get_sample_value(cgbr->tab);
2920
1.00k
2921
1.00k
  
if (1.00k
use_shifted(cgbr)1.00k
) {
2922
191
    struct isl_vec *sample;
2923
191
2924
191
    sample = gbr_get_shifted_sample(cgbr);
2925
191
    if (
!sample || 191
sample->size > 0191
)
2926
42
      return sample;
2927
149
2928
149
    isl_vec_free(sample);
2929
149
  }
2930
1.00k
2931
964
  
if (964
!cgbr->cone964
) {
2932
739
    bset = isl_tab_peek_bset(cgbr->tab);
2933
739
    cgbr->cone = isl_tab_from_recession_cone(bset, 0);
2934
739
    if (!cgbr->cone)
2935
0
      return NULL;
2936
739
    
if (739
isl_tab_track_bset(cgbr->cone,
2937
739
          isl_basic_set_copy(bset)) < 0)
2938
0
      return NULL;
2939
964
  }
2940
964
  
if (964
isl_tab_detect_implicit_equalities(cgbr->cone) < 0964
)
2941
0
    return NULL;
2942
964
2943
964
  
if (964
cgbr->cone->n_dead == cgbr->cone->n_col964
) {
2944
169
    struct isl_vec *sample;
2945
169
    struct isl_tab_undo *snap;
2946
169
2947
169
    if (
cgbr->tab->basis169
) {
2948
48
      if (
cgbr->tab->basis->n_col != 1 + cgbr->tab->n_var48
) {
2949
20
        isl_mat_free(cgbr->tab->basis);
2950
20
        cgbr->tab->basis = NULL;
2951
20
      }
2952
48
      cgbr->tab->n_zero = 0;
2953
48
      cgbr->tab->n_unbounded = 0;
2954
48
    }
2955
169
2956
169
    snap = isl_tab_snap(cgbr->tab);
2957
169
2958
169
    sample = isl_tab_sample(cgbr->tab);
2959
169
2960
169
    if (
!sample || 169
isl_tab_rollback(cgbr->tab, snap) < 0169
) {
2961
0
      isl_vec_free(sample);
2962
0
      return NULL;
2963
0
    }
2964
169
2965
169
    return sample;
2966
169
  }
2967
795
2968
795
  cone = isl_basic_set_dup(isl_tab_peek_bset(cgbr->cone));
2969
795
  cone = drop_constant_terms(cone);
2970
795
  cone = isl_basic_set_update_from_tab(cone, cgbr->cone);
2971
795
  cone = isl_basic_set_underlying_set(cone);
2972
795
  cone = isl_basic_set_gauss(cone, NULL);
2973
795
2974
795
  bset = isl_basic_set_dup(isl_tab_peek_bset(cgbr->tab));
2975
795
  bset = isl_basic_set_update_from_tab(bset, cgbr->tab);
2976
795
  bset = isl_basic_set_underlying_set(bset);
2977
795
  bset = isl_basic_set_gauss(bset, NULL);
2978
795
2979
795
  return isl_basic_set_sample_with_cone(bset, cone);
2980
795
}
2981
2982
static void check_gbr_integer_feasible(struct isl_context_gbr *cgbr)
2983
36.7k
{
2984
36.7k
  struct isl_vec *sample;
2985
36.7k
2986
36.7k
  if (!cgbr->tab)
2987
0
    return;
2988
36.7k
2989
36.7k
  
if (36.7k
cgbr->tab->empty36.7k
)
2990
25.7k
    return;
2991
10.9k
2992
10.9k
  sample = gbr_get_sample(cgbr);
2993
10.9k
  if (!sample)
2994
0
    goto error;
2995
10.9k
2996
10.9k
  
if (10.9k
sample->size == 010.9k
) {
2997
239
    isl_vec_free(sample);
2998
239
    if (isl_tab_mark_empty(cgbr->tab) < 0)
2999
0
      goto error;
3000
239
    return;
3001
239
  }
3002
10.6k
3003
10.6k
  
if (10.6k
isl_tab_add_sample(cgbr->tab, sample) < 010.6k
)
3004
0
    goto error;
3005
10.6k
3006
10.6k
  return;
3007
0
error:
3008
0
  isl_tab_free(cgbr->tab);
3009
0
  cgbr->tab = NULL;
3010
0
}
3011
3012
static struct isl_tab *add_gbr_eq(struct isl_tab *tab, isl_int *eq)
3013
9.12k
{
3014
9.12k
  if (!tab)
3015
0
    return NULL;
3016
9.12k
3017
9.12k
  
if (9.12k
isl_tab_extend_cons(tab, 2) < 09.12k
)
3018
0
    goto error;
3019
9.12k
3020
9.12k
  
if (9.12k
isl_tab_add_eq(tab, eq) < 09.12k
)
3021
0
    goto error;
3022
9.12k
3023
9.12k
  return tab;
3024
0
error:
3025
0
  isl_tab_free(tab);
3026
0
  return NULL;
3027
9.12k
}
3028
3029
/* Add the equality described by "eq" to the context.
3030
 * If "check" is set, then we check if the context is empty after
3031
 * adding the equality.
3032
 * If "update" is set, then we check if the samples are still valid.
3033
 *
3034
 * We do not explicitly add shifted copies of the equality to
3035
 * cgbr->shifted since they would conflict with each other.
3036
 * Instead, we directly mark cgbr->shifted empty.
3037
 */
3038
static void context_gbr_add_eq(struct isl_context *context, isl_int *eq,
3039
    int check, int update)
3040
9.12k
{
3041
9.12k
  struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3042
9.12k
3043
9.12k
  cgbr->tab = add_gbr_eq(cgbr->tab, eq);
3044
9.12k
3045
9.12k
  if (
cgbr->shifted && 9.12k
!cgbr->shifted->empty53
&&
use_shifted(cgbr)1
) {
3046
1
    if (isl_tab_mark_empty(cgbr->shifted) < 0)
3047
0
      goto error;
3048
9.12k
  }
3049
9.12k
3050
9.12k
  
if (9.12k
cgbr->cone && 9.12k
cgbr->cone->n_col != cgbr->cone->n_dead464
) {
3051
442
    if (isl_tab_extend_cons(cgbr->cone, 2) < 0)
3052
0
      goto error;
3053
442
    
if (442
isl_tab_add_eq(cgbr->cone, eq) < 0442
)
3054
0
      goto error;
3055
9.12k
  }
3056
9.12k
3057
9.12k
  
if (9.12k
check9.12k
) {
3058
9.12k
    int v = tab_has_valid_sample(cgbr->tab, eq, 1);
3059
9.12k
    if (v < 0)
3060
0
      goto error;
3061
9.12k
    
if (9.12k
!v9.12k
)
3062
103
      check_gbr_integer_feasible(cgbr);
3063
9.12k
  }
3064
9.12k
  
if (9.12k
update9.12k
)
3065
9.12k
    cgbr->tab = check_samples(cgbr->tab, eq, 1);
3066
9.12k
  return;
3067
0
error:
3068
0
  isl_tab_free(cgbr->tab);
3069
0
  cgbr->tab = NULL;
3070
0
}
3071
3072
static void add_gbr_ineq(struct isl_context_gbr *cgbr, isl_int *ineq)
3073
36.0k
{
3074
36.0k
  if (!cgbr->tab)
3075
0
    return;
3076
36.0k
3077
36.0k
  
if (36.0k
isl_tab_extend_cons(cgbr->tab, 1) < 036.0k
)
3078
0
    goto error;
3079
36.0k
3080
36.0k
  
if (36.0k
isl_tab_add_ineq(cgbr->tab, ineq) < 036.0k
)
3081
0
    goto error;
3082
36.0k
3083
36.0k
  
if (36.0k
cgbr->shifted && 36.0k
!cgbr->shifted->empty461
&&
use_shifted(cgbr)2
) {
3084
2
    int i;
3085
2
    unsigned dim;
3086
2
    dim = isl_basic_map_total_dim(cgbr->tab->bmap);
3087
2
3088
2
    if (isl_tab_extend_cons(cgbr->shifted, 1) < 0)
3089
0
      goto error;
3090
2
3091
10
    
for (i = 0; 2
i < dim10
;
++i8
) {
3092
8
      if (
!8
isl_int_is_neg8
(ineq[1 + i]))
3093
6
        continue;
3094
2
      
isl_int_add2
(ineq[0], ineq[0], ineq[1 + i]);
3095
2
    }
3096
2
3097
2
    if (isl_tab_add_ineq(cgbr->shifted, ineq) < 0)
3098
0
      goto error;
3099
2
3100
10
    
for (i = 0; 2
i < dim10
;
++i8
) {
3101
8
      if (
!8
isl_int_is_neg8
(ineq[1 + i]))
3102
6
        continue;
3103
2
      
isl_int_sub2
(ineq[0], ineq[0], ineq[1 + i]);
3104
2
    }
3105
2
  }
3106
36.0k
3107
36.0k
  
if (36.0k
cgbr->cone && 36.0k
cgbr->cone->n_col != cgbr->cone->n_dead4.30k
) {
3108
3.09k
    if (isl_tab_extend_cons(cgbr->cone, 1) < 0)
3109
0
      goto error;
3110
3.09k
    
if (3.09k
isl_tab_add_ineq(cgbr->cone, ineq) < 03.09k
)
3111
0
      goto error;
3112
36.0k
  }
3113
36.0k
3114
36.0k
  return;
3115
0
error:
3116
0
  isl_tab_free(cgbr->tab);
3117
0
  cgbr->tab = NULL;
3118
0
}
3119
3120
static void context_gbr_add_ineq(struct isl_context *context, isl_int *ineq,
3121
    int check, int update)
3122
25.2k
{
3123
25.2k
  struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3124
25.2k
3125
25.2k
  add_gbr_ineq(cgbr, ineq);
3126
25.2k
  if (!cgbr->tab)
3127
0
    return;
3128
25.2k
3129
25.2k
  
if (25.2k
check25.2k
) {
3130
18.4k
    int v = tab_has_valid_sample(cgbr->tab, ineq, 0);
3131
18.4k
    if (v < 0)
3132
0
      goto error;
3133
18.4k
    
if (18.4k
!v18.4k
)
3134
18.2k
      check_gbr_integer_feasible(cgbr);
3135
18.4k
  }
3136
25.2k
  
if (25.2k
update25.2k
)
3137
5.77k
    cgbr->tab = check_samples(cgbr->tab, ineq, 0);
3138
25.2k
  return;
3139
0
error:
3140
0
  isl_tab_free(cgbr->tab);
3141
0
  cgbr->tab = NULL;
3142
0
}
3143
3144
static isl_stat context_gbr_add_ineq_wrap(void *user, isl_int *ineq)
3145
1.13k
{
3146
1.13k
  struct isl_context *context = (struct isl_context *)user;
3147
1.13k
  context_gbr_add_ineq(context, ineq, 0, 0);
3148
1.13k
  return context->op->is_ok(context) ? 
isl_stat_ok1.13k
:
isl_stat_error0
;
3149
1.13k
}
3150
3151
static enum isl_tab_row_sign context_gbr_ineq_sign(struct isl_context *context,
3152
      isl_int *ineq, int strict)
3153
11.3k
{
3154
11.3k
  struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3155
11.3k
  return tab_ineq_sign(cgbr->tab, ineq, strict);
3156
11.3k
}
3157
3158
/* Check whether "ineq" can be added to the tableau without rendering
3159
 * it infeasible.
3160
 */
3161
static int context_gbr_test_ineq(struct isl_context *context, isl_int *ineq)
3162
10.7k
{
3163
10.7k
  struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3164
10.7k
  struct isl_tab_undo *snap;
3165
10.7k
  struct isl_tab_undo *shifted_snap = NULL;
3166
10.7k
  struct isl_tab_undo *cone_snap = NULL;
3167
10.7k
  int feasible;
3168
10.7k
3169
10.7k
  if (!cgbr->tab)
3170
0
    return -1;
3171
10.7k
3172
10.7k
  
if (10.7k
isl_tab_extend_cons(cgbr->tab, 1) < 010.7k
)
3173
0
    return -1;
3174
10.7k
3175
10.7k
  snap = isl_tab_snap(cgbr->tab);
3176
10.7k
  if (cgbr->shifted)
3177
311
    shifted_snap = isl_tab_snap(cgbr->shifted);
3178
10.7k
  if (cgbr->cone)
3179
2.16k
    cone_snap = isl_tab_snap(cgbr->cone);
3180
10.7k
  add_gbr_ineq(cgbr, ineq);
3181
10.7k
  check_gbr_integer_feasible(cgbr);
3182
10.7k
  if (!cgbr->tab)
3183
0
    return -1;
3184
10.7k
  feasible = !cgbr->tab->empty;
3185
10.7k
  if (isl_tab_rollback(cgbr->tab, snap) < 0)
3186
0
    return -1;
3187
10.7k
  
if (10.7k
shifted_snap10.7k
) {
3188
311
    if (isl_tab_rollback(cgbr->shifted, shifted_snap))
3189
0
      return -1;
3190
10.4k
  } else 
if (10.4k
cgbr->shifted10.4k
) {
3191
104
    isl_tab_free(cgbr->shifted);
3192
104
    cgbr->shifted = NULL;
3193
104
  }
3194
10.7k
  
if (10.7k
cone_snap10.7k
) {
3195
2.16k
    if (isl_tab_rollback(cgbr->cone, cone_snap))
3196
0
      return -1;
3197
8.62k
  } else 
if (8.62k
cgbr->cone8.62k
) {
3198
280
    isl_tab_free(cgbr->cone);
3199
280
    cgbr->cone = NULL;
3200
280
  }
3201
10.7k
3202
10.7k
  return feasible;
3203
10.7k
}
3204
3205
/* Return the column of the last of the variables associated to
3206
 * a column that has a non-zero coefficient.
3207
 * This function is called in a context where only coefficients
3208
 * of parameters or divs can be non-zero.
3209
 */
3210
static int last_non_zero_var_col(struct isl_tab *tab, isl_int *p)
3211
163
{
3212
163
  int i;
3213
163
  int col;
3214
163
3215
163
  if (tab->n_var == 0)
3216
0
    return -1;
3217
163
3218
449
  
for (i = tab->n_var - 1; 163
i >= 0449
;
--i286
) {
3219
449
    if (
i >= tab->n_param && 449
i < tab->n_var - tab->n_div432
)
3220
31
      continue;
3221
418
    
if (418
tab->var[i].is_row418
)
3222
72
      continue;
3223
346
    col = tab->var[i].index;
3224
346
    if (
!346
isl_int_is_zero346
(p[col]))
3225
163
      return col;
3226
449
  }
3227
163
3228
0
  return -1;
3229
163
}
3230
3231
/* Look through all the recently added equalities in the context
3232
 * to see if we can propagate any of them to the main tableau.
3233
 *
3234
 * The newly added equalities in the context are encoded as pairs
3235
 * of inequalities starting at inequality "first".
3236
 *
3237
 * We tentatively add each of these equalities to the main tableau
3238
 * and if this happens to result in a row with a final coefficient
3239
 * that is one or negative one, we use it to kill a column
3240
 * in the main tableau.  Otherwise, we discard the tentatively
3241
 * added row.
3242
 * This tentative addition of equality constraints turns
3243
 * on the undo facility of the tableau.  Turn it off again
3244
 * at the end, assuming it was turned off to begin with.
3245
 *
3246
 * Return 0 on success and -1 on failure.
3247
 */
3248
static int propagate_equalities(struct isl_context_gbr *cgbr,
3249
  struct isl_tab *tab, unsigned first)
3250
102
{
3251
102
  int i;
3252
102
  struct isl_vec *eq = NULL;
3253
102
  isl_bool needs_undo;
3254
102
3255
102
  needs_undo = isl_tab_need_undo(tab);
3256
102
  if (needs_undo < 0)
3257
0
    goto error;
3258
102
  eq = isl_vec_alloc(tab->mat->ctx, 1 + tab->n_var);
3259
102
  if (!eq)
3260
0
    goto error;
3261
102
3262
102
  
if (102
isl_tab_extend_cons(tab, (cgbr->tab->bmap->n_ineq - first)/2) < 0102
)
3263
0
    goto error;
3264
102
3265
102
  isl_seq_clr(eq->el + 1 + tab->n_param,
3266
102
        tab->n_var - tab->n_param - tab->n_div);
3267
265
  for (i = first; 
i < cgbr->tab->bmap->n_ineq265
;
i += 2163
) {
3268
163
    int j;
3269
163
    int r;
3270
163
    struct isl_tab_undo *snap;
3271
163
    snap = isl_tab_snap(tab);
3272
163
3273
163
    isl_seq_cpy(eq->el, cgbr->tab->bmap->ineq[i], 1 + tab->n_param);
3274
163
    isl_seq_cpy(eq->el + 1 + tab->n_var - tab->n_div,
3275
163
          cgbr->tab->bmap->ineq[i] + 1 + tab->n_param,
3276
163
          tab->n_div);
3277
163
3278
163
    r = isl_tab_add_row(tab, eq->el);
3279
163
    if (r < 0)
3280
0
      goto error;
3281
163
    r = tab->con[r].index;
3282
163
    j = last_non_zero_var_col(tab, tab->mat->row[r] + 2 + tab->M);
3283
163
    if (
j < 0 || 163
j < tab->n_dead163
||
3284
163
        
!163
isl_int_is_one163
(tab->mat->row[r][0]) ||
3285
163
        
(!163
isl_int_is_one163
(tab->mat->row[r][2 + tab->M + j]) &&
3286
163
         
!100
isl_int_is_negone100
(tab->mat->row[r][2 + tab->M + j]))) {
3287
53
      if (isl_tab_rollback(tab, snap) < 0)
3288
0
        goto error;
3289
53
      continue;
3290
53
    }
3291
110
    
if (110
isl_tab_pivot(tab, r, j) < 0110
)
3292
0
      goto error;
3293
110
    
if (110
isl_tab_kill_col(tab, j) < 0110
)
3294
0
      goto error;
3295
110
3296
110
    
if (110
restore_lexmin(tab) < 0110
)
3297
0
      goto error;
3298
163
  }
3299
102
3300
102
  
if (102
!needs_undo102
)
3301
102
    isl_tab_clear_undo(tab);
3302
102
  isl_vec_free(eq);
3303
102
3304
102
  return 0;
3305
0
error:
3306
0
  isl_vec_free(eq);
3307
0
  isl_tab_free(cgbr->tab);
3308
0
  cgbr->tab = NULL;
3309
0
  return -1;
3310
102
}
3311
3312
static int context_gbr_detect_equalities(struct isl_context *context,
3313
  struct isl_tab *tab)
3314
502
{
3315
502
  struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3316
502
  unsigned n_ineq;
3317
502
3318
502
  if (
!cgbr->cone502
) {
3319
226
    struct isl_basic_set *bset = isl_tab_peek_bset(cgbr->tab);
3320
226
    cgbr->cone = isl_tab_from_recession_cone(bset, 0);
3321
226
    if (!cgbr->cone)
3322
0
      goto error;
3323
226
    
if (226
isl_tab_track_bset(cgbr->cone,
3324
226
          isl_basic_set_copy(bset)) < 0)
3325
0
      goto error;
3326
502
  }
3327
502
  
if (502
isl_tab_detect_implicit_equalities(cgbr->cone) < 0502
)
3328
0
    goto error;
3329
502
3330
502
  n_ineq = cgbr->tab->bmap->n_ineq;
3331
502
  cgbr->tab = isl_tab_detect_equalities(cgbr->tab, cgbr->cone);
3332
502
  if (!cgbr->tab)
3333
0
    return -1;
3334
502
  
if (502
cgbr->tab->bmap->n_ineq > n_ineq &&
3335
102
      propagate_equalities(cgbr, tab, n_ineq) < 0)
3336
0
    return -1;
3337
502
3338
502
  return 0;
3339
0
error:
3340
0
  isl_tab_free(cgbr->tab);
3341
0
  cgbr->tab = NULL;
3342
0
  return -1;
3343
502
}
3344
3345
static int context_gbr_get_div(struct isl_context *context, struct isl_tab *tab,
3346
    struct isl_vec *div)
3347
592
{
3348
592
  return get_div(tab, context, div);
3349
592
}
3350
3351
static isl_bool context_gbr_insert_div(struct isl_context *context, int pos,
3352
  __isl_keep isl_vec *div)
3353
565
{
3354
565
  struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3355
565
  if (
cgbr->cone565
) {
3356
283
    int r, n_div, o_div;
3357
283
3358
283
    n_div = isl_basic_map_dim(cgbr->cone->bmap, isl_dim_div);
3359
283
    o_div = cgbr->cone->n_var - n_div;
3360
283
3361
283
    if (isl_tab_extend_cons(cgbr->cone, 3) < 0)
3362
0
      return isl_bool_error;
3363
283
    
if (283
isl_tab_extend_vars(cgbr->cone, 1) < 0283
)
3364
0
      return isl_bool_error;
3365
283
    
if (283
(r = isl_tab_insert_var(cgbr->cone, pos)) <0283
)
3366
0
      return isl_bool_error;
3367
283
3368
283
    cgbr->cone->bmap = isl_basic_map_insert_div(cgbr->cone->bmap,
3369
283
                r - o_div, div);
3370
283
    if (!cgbr->cone->bmap)
3371
0
      return isl_bool_error;
3372
283
    
if (283
isl_tab_push_var(cgbr->cone, isl_tab_undo_bmap_div,
3373
283
            &cgbr->cone->var[r]) < 0)
3374
0
      return isl_bool_error;
3375
565
  }
3376
565
  return context_tab_insert_div(cgbr->tab, pos, div,
3377
565
          context_gbr_add_ineq_wrap, context);
3378
565
}
3379
3380
static int context_gbr_best_split(struct isl_context *context,
3381
    struct isl_tab *tab)
3382
331
{
3383
331
  struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3384
331
  struct isl_tab_undo *snap;
3385
331
  int r;
3386
331
3387
331
  snap = isl_tab_snap(cgbr->tab);
3388
331
  r = best_split(tab, cgbr->tab);
3389
331
3390
331
  if (
r >= 0 && 331
isl_tab_rollback(cgbr->tab, snap) < 0331
)
3391
0
    return -1;
3392
331
3393
331
  return r;
3394
331
}
3395
3396
static int context_gbr_is_empty(struct isl_context *context)
3397
43.8k
{
3398
43.8k
  struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3399
43.8k
  if (!cgbr->tab)
3400
0
    return -1;
3401
43.8k
  return cgbr->tab->empty;
3402
43.8k
}
3403
3404
struct isl_gbr_tab_undo {
3405
  struct isl_tab_undo *tab_snap;
3406
  struct isl_tab_undo *shifted_snap;
3407
  struct isl_tab_undo *cone_snap;
3408
};
3409
3410
static void *context_gbr_save(struct isl_context *context)
3411
27.7k
{
3412
27.7k
  struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3413
27.7k
  struct isl_gbr_tab_undo *snap;
3414
27.7k
3415
27.7k
  if (!cgbr->tab)
3416
0
    return NULL;
3417
27.7k
3418
27.7k
  
snap = 27.7k
isl_alloc_type27.7k
(cgbr->tab->mat->ctx, struct isl_gbr_tab_undo);
3419
27.7k
  if (!snap)
3420
0
    return NULL;
3421
27.7k
3422
27.7k
  snap->tab_snap = isl_tab_snap(cgbr->tab);
3423
27.7k
  if (isl_tab_save_samples(cgbr->tab) < 0)
3424
0
    goto error;
3425
27.7k
3426
27.7k
  
if (27.7k
cgbr->shifted27.7k
)
3427
163
    snap->shifted_snap = isl_tab_snap(cgbr->shifted);
3428
27.7k
  else
3429
27.5k
    snap->shifted_snap = NULL;
3430
27.7k
3431
27.7k
  if (cgbr->cone)
3432
1.45k
    snap->cone_snap = isl_tab_snap(cgbr->cone);
3433
27.7k
  else
3434
26.2k
    snap->cone_snap = NULL;
3435
27.7k
3436
27.7k
  return snap;
3437
0
error:
3438
0
  free(snap);
3439
0
  return NULL;
3440
27.7k
}
3441
3442
static void context_gbr_restore(struct isl_context *context, void *save)
3443
23.0k
{
3444
23.0k
  struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3445
23.0k
  struct isl_gbr_tab_undo *snap = (struct isl_gbr_tab_undo *)save;
3446
23.0k
  if (!snap)
3447
0
    goto error;
3448
23.0k
  
if (23.0k
isl_tab_rollback(cgbr->tab, snap->tab_snap) < 023.0k
)
3449
0
    goto error;
3450
23.0k
3451
23.0k
  
if (23.0k
snap->shifted_snap23.0k
) {
3452
115
    if (isl_tab_rollback(cgbr->shifted, snap->shifted_snap) < 0)
3453
0
      goto error;
3454
22.9k
  } else 
if (22.9k
cgbr->shifted22.9k
) {
3455
0
    isl_tab_free(cgbr->shifted);
3456
0
    cgbr->shifted = NULL;
3457
0
  }
3458
23.0k
3459
23.0k
  
if (23.0k
snap->cone_snap23.0k
) {
3460
1.38k
    if (isl_tab_rollback(cgbr->cone, snap->cone_snap) < 0)
3461
0
      goto error;
3462
21.6k
  } else 
if (21.6k
cgbr->cone21.6k
) {
3463
250
    isl_tab_free(cgbr->cone);
3464
250
    cgbr->cone = NULL;
3465
250
  }
3466
23.0k
3467
23.0k
  free(snap);
3468
23.0k
3469
23.0k
  return;
3470
0
error:
3471
0
  free(snap);
3472
0
  isl_tab_free(cgbr->tab);
3473
0
  cgbr->tab = NULL;
3474
0
}
3475
3476
static void context_gbr_discard(void *save)
3477
4.65k
{
3478
4.65k
  struct isl_gbr_tab_undo *snap = (struct isl_gbr_tab_undo *)save;
3479
4.65k
  free(snap);
3480
4.65k
}
3481
3482
static int context_gbr_is_ok(struct isl_context *context)
3483
1.78k
{
3484
1.78k
  struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3485
1.78k
  return !!cgbr->tab;
3486
1.78k
}
3487
3488
static void context_gbr_invalidate(struct isl_context *context)
3489
0
{
3490
0
  struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3491
0
  isl_tab_free(cgbr->tab);
3492
0
  cgbr->tab = NULL;
3493
0
}
3494
3495
static __isl_null struct isl_context *context_gbr_free(
3496
  struct isl_context *context)
3497
7.57k
{
3498
7.57k
  struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3499
7.57k
  isl_tab_free(cgbr->tab);
3500
7.57k
  isl_tab_free(cgbr->shifted);
3501
7.57k
  isl_tab_free(cgbr->cone);
3502
7.57k
  free(cgbr);
3503
7.57k
3504
7.57k
  return NULL;
3505
7.57k
}
3506
3507
struct isl_context_op isl_context_gbr_op = {
3508
  context_gbr_detect_nonnegative_parameters,
3509
  context_gbr_peek_basic_set,
3510
  context_gbr_peek_tab,
3511
  context_gbr_add_eq,
3512
  context_gbr_add_ineq,
3513
  context_gbr_ineq_sign,
3514
  context_gbr_test_ineq,
3515
  context_gbr_get_div,
3516
  context_gbr_insert_div,
3517
  context_gbr_detect_equalities,
3518
  context_gbr_best_split,
3519
  context_gbr_is_empty,
3520
  context_gbr_is_ok,
3521
  context_gbr_save,
3522
  context_gbr_restore,
3523
  context_gbr_discard,
3524
  context_gbr_invalidate,
3525
  context_gbr_free,
3526
};
3527
3528
static struct isl_context *isl_context_gbr_alloc(__isl_keep isl_basic_set *dom)
3529
7.57k
{
3530
7.57k
  struct isl_context_gbr *cgbr;
3531
7.57k
3532
7.57k
  if (!dom)
3533
0
    return NULL;
3534
7.57k
3535
7.57k
  
cgbr = 7.57k
isl_calloc_type7.57k
(dom->ctx, struct isl_context_gbr);
3536
7.57k
  if (!cgbr)
3537
0
    return NULL;
3538
7.57k
3539
7.57k
  cgbr->context.op = &isl_context_gbr_op;
3540
7.57k
3541
7.57k
  cgbr->shifted = NULL;
3542
7.57k
  cgbr->cone = NULL;
3543
7.57k
  cgbr->tab = isl_tab_from_basic_set(dom, 1);
3544
7.57k
  cgbr->tab = isl_tab_init_samples(cgbr->tab);
3545
7.57k
  if (!cgbr->tab)
3546
0
    goto error;
3547
7.57k
  check_gbr_integer_feasible(cgbr);
3548
7.57k
3549
7.57k
  return &cgbr->context;
3550
0
error:
3551
0
  cgbr->context.op->free(&cgbr->context);
3552
0
  return NULL;
3553
7.57k
}
3554
3555
/* Allocate a context corresponding to "dom".
3556
 * The representation specific fields are initialized by
3557
 * isl_context_lex_alloc or isl_context_gbr_alloc.
3558
 * The shared "n_unknown" field is initialized to the number
3559
 * of final unknown integer divisions in "dom".
3560
 */
3561
static struct isl_context *isl_context_alloc(__isl_keep isl_basic_set *dom)
3562
7.57k
{
3563
7.57k
  struct isl_context *context;
3564
7.57k
  int first;
3565
7.57k
3566
7.57k
  if (!dom)
3567
0
    return NULL;
3568
7.57k
3569
7.57k
  
if (7.57k
dom->ctx->opt->context == 7.57k
ISL_CONTEXT_LEXMIN7.57k
)
3570
0
    context = isl_context_lex_alloc(dom);
3571
7.57k
  else
3572
7.57k
    context = isl_context_gbr_alloc(dom);
3573
7.57k
3574
7.57k
  if (!context)
3575
0
    return NULL;
3576
7.57k
3577
7.57k
  first = isl_basic_set_first_unknown_div(dom);
3578
7.57k
  if (first < 0)
3579
0
    return context->op->free(context);
3580
7.57k
  context->n_unknown = isl_basic_set_dim(dom, isl_dim_div) - first;
3581
7.57k
3582
7.57k
  return context;
3583
7.57k
}
3584
3585
/* Initialize some common fields of "sol", which keeps track
3586
 * of the solution of an optimization problem on "bmap" over
3587
 * the domain "dom".
3588
 * If "max" is set, then a maximization problem is being solved, rather than
3589
 * a minimization problem, which means that the variables in the
3590
 * tableau have value "M - x" rather than "M + x".
3591
 */
3592
static isl_stat sol_init(struct isl_sol *sol, __isl_keep isl_basic_map *bmap,
3593
  __isl_keep isl_basic_set *dom, int max)
3594
7.57k
{
3595
7.57k
  sol->rational = ISL_F_ISSET(bmap, ISL_BASIC_MAP_RATIONAL);
3596
7.57k
  sol->dec_level.callback.run = &sol_dec_level_wrap;
3597
7.57k
  sol->dec_level.sol = sol;
3598
7.57k
  sol->max = max;
3599
7.57k
  sol->n_out = isl_basic_map_dim(bmap, isl_dim_out);
3600
7.57k
  sol->space = isl_basic_map_get_space(bmap);
3601
7.57k
3602
7.57k
  sol->context = isl_context_alloc(dom);
3603
7.57k
  if (
!sol->space || 7.57k
!sol->context7.57k
)
3604
0
    return isl_stat_error;
3605
7.57k
3606
7.57k
  return isl_stat_ok;
3607
7.57k
}
3608
3609
/* Construct an isl_sol_map structure for accumulating the solution.
3610
 * If track_empty is set, then we also keep track of the parts
3611
 * of the context where there is no solution.
3612
 * If max is set, then we are solving a maximization, rather than
3613
 * a minimization problem, which means that the variables in the
3614
 * tableau have value "M - x" rather than "M + x".
3615
 */
3616
static struct isl_sol *sol_map_init(__isl_keep isl_basic_map *bmap,
3617
  __isl_take isl_basic_set *dom, int track_empty, int max)
3618
3.40k
{
3619
3.40k
  struct isl_sol_map *sol_map = NULL;
3620
3.40k
  isl_space *space;
3621
3.40k
3622
3.40k
  if (!bmap)
3623
0
    goto error;
3624
3.40k
3625
3.40k
  
sol_map = 3.40k
isl_calloc_type3.40k
(bmap->ctx, struct isl_sol_map);
3626
3.40k
  if (!sol_map)
3627
0
    goto error;
3628
3.40k
3629
3.40k
  sol_map->sol.free = &sol_map_free;
3630
3.40k
  if (sol_init(&sol_map->sol, bmap, dom, max) < 0)
3631
0
    goto error;
3632
3.40k
  sol_map->sol.add = &sol_map_add_wrap;
3633
2.60k
  sol_map->sol.add_empty = track_empty ? &sol_map_add_empty_wrap : NULL;
3634
3.40k
  space = isl_space_copy(sol_map->sol.space);
3635
3.40k
  sol_map->map = isl_map_alloc_space(space, 1, ISL_MAP_DISJOINT);
3636
3.40k
  if (!sol_map->map)
3637
0
    goto error;
3638
3.40k
3639
3.40k
  
if (3.40k
track_empty3.40k
) {
3640
2.60k
    sol_map->empty = isl_set_alloc_space(isl_basic_set_get_space(dom),
3641
2.60k
              1, ISL_SET_DISJOINT);
3642
2.60k
    if (!sol_map->empty)
3643
0
      goto error;
3644
3.40k
  }
3645
3.40k
3646
3.40k
  isl_basic_set_free(dom);
3647
3.40k
  return &sol_map->sol;
3648
0
error:
3649
0
  isl_basic_set_free(dom);
3650
0
  sol_free(&sol_map->sol);
3651
0
  return NULL;
3652
3.40k
}
3653
3654
/* Check whether all coefficients of (non-parameter) variables
3655
 * are non-positive, meaning that no pivots can be performed on the row.
3656
 */
3657
static int is_critical(struct isl_tab *tab, int row)
3658
11.3k
{
3659
11.3k
  int j;
3660
11.3k
  unsigned off = 2 + tab->M;
3661
11.3k
3662
48.7k
  for (j = tab->n_dead; 
j < tab->n_col48.7k
;
++j37.4k
) {
3663
38.8k
    if (tab->col_var[j] >= 0 &&
3664
27.0k
        (tab->col_var[j] < tab->n_param  ||
3665
1.45k
        tab->col_var[j] >= tab->n_var - tab->n_div))
3666
27.0k
      continue;
3667
11.7k
3668
11.7k
    
if (11.7k
isl_int_is_pos11.7k
(tab->mat->row[row][off + j]))
3669
1.36k
      return 0;
3670
38.8k
  }
3671
11.3k
3672
9.94k
  return 1;
3673
11.3k
}
3674
3675
/* Check whether the inequality represented by vec is strict over the integers,
3676
 * i.e., there are no integer values satisfying the constraint with
3677
 * equality.  This happens if the gcd of the coefficients is not a divisor
3678
 * of the constant term.  If so, scale the constraint down by the gcd
3679
 * of the coefficients.
3680
 */
3681
static int is_strict(struct isl_vec *vec)
3682
14.1k
{
3683
14.1k
  isl_int gcd;
3684
14.1k
  int strict = 0;
3685
14.1k
3686
14.1k
  isl_int_init(gcd);
3687
14.1k
  isl_seq_gcd(vec->el + 1, vec->size - 1, &gcd);
3688
14.1k
  if (
!14.1k
isl_int_is_one14.1k
(gcd)) {
3689
1.09k
    strict = !isl_int_is_divisible_by(vec->el[0], gcd);
3690
1.09k
    isl_int_fdiv_q(vec->el[0], vec->el[0], gcd);
3691
1.09k
    isl_seq_scale_down(vec->el + 1, vec->el + 1, gcd, vec->size-1);
3692
1.09k
  }
3693
14.1k
  isl_int_clear(gcd);
3694
14.1k
3695
14.1k
  return strict;
3696
14.1k
}
3697
3698
/* Determine the sign of the given row of the main tableau.
3699
 * The result is one of
3700
 *  isl_tab_row_pos: always non-negative; no pivot needed
3701
 *  isl_tab_row_neg: always non-positive; pivot
3702
 *  isl_tab_row_any: can be both positive and negative; split
3703
 *
3704
 * We first handle some simple cases
3705
 *  - the row sign may be known already
3706
 *  - the row may be obviously non-negative
3707
 *  - the parametric constant may be equal to that of another row
3708
 *    for which we know the sign.  This sign will be either "pos" or
3709
 *    "any".  If it had been "neg" then we would have pivoted before.
3710
 *
3711
 * If none of these cases hold, we check the value of the row for each
3712
 * of the currently active samples.  Based on the signs of these values
3713
 * we make an initial determination of the sign of the row.
3714
 *
3715
 *  all zero      ->  unk(nown)
3716
 *  all non-negative    ->  pos
3717
 *  all non-positive    ->  neg
3718
 *  both negative and positive  ->  all
3719
 *
3720
 * If we end up with "all", we are done.
3721
 * Otherwise, we perform a check for positive and/or negative
3722
 * values as follows.
3723
 *
3724
 *  samples        neg         unk         pos
3725
 *  <0 ?          Y        N      Y        N
3726
 *              pos    any      pos
3727
 *  >0 ?       Y      N  Y     N
3728
 *        any    neg  any   neg
3729
 *
3730
 * There is no special sign for "zero", because we can usually treat zero
3731
 * as either non-negative or non-positive, whatever works out best.
3732
 * However, if the row is "critical", meaning that pivoting is impossible
3733
 * then we don't want to limp zero with the non-positive case, because
3734
 * then we we would lose the solution for those values of the parameters
3735
 * where the value of the row is zero.  Instead, we treat 0 as non-negative
3736
 * ensuring a split if the row can attain both zero and negative values.
3737
 * The same happens when the original constraint was one that could not
3738
 * be satisfied with equality by any integer values of the parameters.
3739
 * In this case, we normalize the constraint, but then a value of zero
3740
 * for the normalized constraint is actually a positive value for the
3741
 * original constraint, so again we need to treat zero as non-negative.
3742
 * In both these cases, we have the following decision tree instead:
3743
 *
3744
 *  all non-negative    ->  pos
3745
 *  all negative      ->  neg
3746
 *  both negative and non-negative  ->  all
3747
 *
3748
 *  samples        neg                     pos
3749
 *  <0 ?                        Y        N
3750
 *                     any      pos
3751
 *  >=0 ?      Y      N
3752
 *        any    neg
3753
 */
3754
static enum isl_tab_row_sign row_sign(struct isl_tab *tab,
3755
  struct isl_sol *sol, int row)
3756
61.4k
{
3757
61.4k
  struct isl_vec *ineq = NULL;
3758
61.4k
  enum isl_tab_row_sign res = isl_tab_row_unknown;
3759
61.4k
  int critical;
3760
61.4k
  int strict;
3761
61.4k
  int row2;
3762
61.4k
3763
61.4k
  if (tab->row_sign[row] != isl_tab_row_unknown)
3764
31.4k
    return tab->row_sign[row];
3765
30.0k
  
if (30.0k
is_obviously_nonneg(tab, row)30.0k
)
3766
18.6k
    return isl_tab_row_pos;
3767
104k
  
for (row2 = tab->n_redundant; 11.3k
row2 < tab->n_row104k
;
++row292.9k
) {
3768
93.0k
    if (tab->row_sign[row2] == isl_tab_row_unknown)
3769
30.3k
      continue;
3770
62.6k
    
if (62.6k
identical_parameter_line(tab, row, row2)62.6k
)
3771
43
      return tab->row_sign[row2];
3772
93.0k
  }
3773
11.3k
3774
11.3k
  critical = is_critical(tab, row);
3775
11.3k
3776
11.3k
  ineq = get_row_parameter_ineq(tab, row);
3777
11.3k
  if (!ineq)
3778
0
    goto error;
3779
11.3k
3780
11.3k
  strict = is_strict(ineq);
3781
11.3k
3782
11.3k
  res = sol->context->op->ineq_sign(sol->context, ineq->el,
3783
1.36k
            critical || strict);
3784
11.3k
3785
11.3k
  if (
res == isl_tab_row_unknown || 11.3k
res == isl_tab_row_pos11.2k
) {
3786
8.64k
    /* test for negative values */
3787
8.64k
    int feasible;
3788
8.64k
    isl_seq_neg(ineq->el, ineq->el, ineq->size);
3789
8.64k
    isl_int_sub_ui(ineq->el[0], ineq->el[0], 1);
3790
8.64k
3791
8.64k
    feasible = sol->context->op->test_ineq(sol->context, ineq->el);
3792
8.64k
    if (feasible < 0)
3793
0
      goto error;
3794
8.64k
    
if (8.64k
!feasible8.64k
)
3795
7.99k
      res = isl_tab_row_pos;
3796
8.64k
    else
3797
652
      
res = (res == isl_tab_row_unknown) ? 652
isl_tab_row_neg15
3798
637
                 : isl_tab_row_any;
3799
8.64k
    if (
res == isl_tab_row_neg8.64k
) {
3800
15
      isl_seq_neg(ineq->el, ineq->el, ineq->size);
3801
15
      isl_int_sub_ui(ineq->el[0], ineq->el[0], 1);
3802
15
    }
3803
8.64k
  }
3804
11.3k
3805
11.3k
  
if (11.3k
res == isl_tab_row_neg11.3k
) {
3806
2.14k
    /* test for positive values */
3807
2.14k
    int feasible;
3808
2.14k
    if (
!critical && 2.14k
!strict81
)
3809
81
      isl_int_sub_ui(ineq->el[0], ineq->el[0], 1);
3810
2.14k
3811
2.14k
    feasible = sol->context->op->test_ineq(sol->context, ineq->el);
3812
2.14k
    if (feasible < 0)
3813
0
      goto error;
3814
2.14k
    
if (2.14k
feasible2.14k
)
3815
2.12k
      res = isl_tab_row_any;
3816
2.14k
  }
3817
11.3k
3818
11.3k
  isl_vec_free(ineq);
3819
11.3k
  return res;
3820
0
error:
3821
0
  isl_vec_free(ineq);
3822
0
  return isl_tab_row_unknown;
3823
61.4k
}
3824
3825
static void find_solutions(struct isl_sol *sol, struct isl_tab *tab);
3826
3827
/* Find solutions for values of the parameters that satisfy the given
3828
 * inequality.
3829
 *
3830
 * We currently take a snapshot of the context tableau that is reset
3831
 * when we return from this function, while we make a copy of the main
3832
 * tableau, leaving the original main tableau untouched.
3833
 * These are fairly arbitrary choices.  Making a copy also of the context
3834
 * tableau would obviate the need to undo any changes made to it later,
3835
 * while taking a snapshot of the main tableau could reduce memory usage.
3836
 * If we were to switch to taking a snapshot of the main tableau,
3837
 * we would have to keep in mind that we need to save the row signs
3838
 * and that we need to do this before saving the current basis
3839
 * such that the basis has been restore before we restore the row signs.
3840
 */
3841
static void find_in_pos(struct isl_sol *sol, struct isl_tab *tab, isl_int *ineq)
3842
2.84k
{
3843
2.84k
  void *saved;
3844
2.84k
3845
2.84k
  if (!sol->context)
3846
0
    goto error;
3847
2.84k
  saved = sol->context->op->save(sol->context);
3848
2.84k
3849
2.84k
  tab = isl_tab_dup(tab);
3850
2.84k
  if (!tab)
3851
0
    goto error;
3852
2.84k
3853
2.84k
  sol->context->op->add_ineq(sol->context, ineq, 0, 1);
3854
2.84k
3855
2.84k
  find_solutions(sol, tab);
3856
2.84k
3857
2.84k
  if (!sol->error)
3858
2.84k
    sol->context->op->restore(sol->context, saved);
3859
2.84k
  else
3860
0
    sol->context->op->discard(saved);
3861
2.84k
  return;
3862
0
error:
3863
0
  sol->error = 1;
3864
0
}
3865
3866
/* Record the absence of solutions for those values of the parameters
3867
 * that do not satisfy the given inequality with equality.
3868
 */
3869
static void no_sol_in_strict(struct isl_sol *sol,
3870
  struct isl_tab *tab, struct isl_vec *ineq)
3871
18.3k
{
3872
18.3k
  int empty;
3873
18.3k
  void *saved;
3874
18.3k
3875
18.3k
  if (
!sol->context || 18.3k
sol->error18.3k
)
3876
0
    goto error;
3877
18.3k
  saved = sol->context->op->save(sol->context);
3878
18.3k
3879
18.3k
  isl_int_sub_ui(ineq->el[0], ineq->el[0], 1);
3880
18.3k
3881
18.3k
  sol->context->op->add_ineq(sol->context, ineq->el, 1, 0);
3882
18.3k
  if (!sol->context)
3883
0
    goto error;
3884
18.3k
3885
18.3k
  empty = tab->empty;
3886
18.3k
  tab->empty = 1;
3887
18.3k
  sol_add(sol, tab);
3888
18.3k
  tab->empty = empty;
3889
18.3k
3890
18.3k
  isl_int_add_ui(ineq->el[0], ineq->el[0], 1);
3891
18.3k
3892
18.3k
  sol->context->op->restore(sol->context, saved);
3893
18.3k
  return;
3894
0
error:
3895
0
  sol->error = 1;
3896
0
}
3897
3898
/* Reset all row variables that are marked to have a sign that may
3899
 * be both positive and negative to have an unknown sign.
3900
 */
3901
static void reset_any_to_unknown(struct isl_tab *tab)
3902
2.84k
{
3903
2.84k
  int row;
3904
2.84k
3905
22.2k
  for (row = tab->n_redundant; 
row < tab->n_row22.2k
;
++row19.4k
) {
3906
19.4k
    if (!isl_tab_var_from_row(tab, row)->is_nonneg)
3907
17
      continue;
3908
19.4k
    
if (19.4k
tab->row_sign[row] == isl_tab_row_any19.4k
)
3909
3.30k
      tab->row_sign[row] = isl_tab_row_unknown;
3910
19.4k
  }
3911
2.84k
}
3912
3913
/* Compute the lexicographic minimum of the set represented by the main
3914
 * tableau "tab" within the context "sol->context_tab".
3915
 * On entry the sample value of the main tableau is lexicographically
3916
 * less than or equal to this lexicographic minimum.
3917
 * Pivots are performed until a feasible point is found, which is then
3918
 * necessarily equal to the minimum, or until the tableau is found to
3919
 * be infeasible.  Some pivots may need to be performed for only some
3920
 * feasible values of the context tableau.  If so, the context tableau
3921
 * is split into a part where the pivot is needed and a part where it is not.
3922
 *
3923
 * Whenever we enter the main loop, the main tableau is such that no
3924
 * "obvious" pivots need to be performed on it, where "obvious" means
3925
 * that the given row can be seen to be negative without looking at
3926
 * the context tableau.  In particular, for non-parametric problems,
3927
 * no pivots need to be performed on the main tableau.
3928
 * The caller of find_solutions is responsible for making this property
3929
 * hold prior to the first iteration of the loop, while restore_lexmin
3930
 * is called before every other iteration.
3931
 *
3932
 * Inside the main loop, we first examine the signs of the rows of
3933
 * the main tableau within the context of the context tableau.
3934
 * If we find a row that is always non-positive for all values of
3935
 * the parameters satisfying the context tableau and negative for at
3936
 * least one value of the parameters, we perform the appropriate pivot
3937
 * and start over.  An exception is the case where no pivot can be
3938
 * performed on the row.  In this case, we require that the sign of
3939
 * the row is negative for all values of the parameters (rather than just
3940
 * non-positive).  This special case is handled inside row_sign, which
3941
 * will say that the row can have any sign if it determines that it can
3942
 * attain both negative and zero values.
3943
 *
3944
 * If we can't find a row that always requires a pivot, but we can find
3945
 * one or more rows that require a pivot for some values of the parameters
3946
 * (i.e., the row can attain both positive and negative signs), then we split
3947
 * the context tableau into two parts, one where we force the sign to be
3948
 * non-negative and one where we force is to be negative.
3949
 * The non-negative part is handled by a recursive call (through find_in_pos).
3950
 * Upon returning from this call, we continue with the negative part and
3951
 * perform the required pivot.
3952
 *
3953
 * If no such rows can be found, all rows are non-negative and we have
3954
 * found a (rational) feasible point.  If we only wanted a rational point
3955
 * then we are done.
3956
 * Otherwise, we check if all values of the sample point of the tableau
3957
 * are integral for the variables.  If so, we have found the minimal
3958
 * integral point and we are done.
3959
 * If the sample point is not integral, then we need to make a distinction
3960
 * based on whether the constant term is non-integral or the coefficients
3961
 * of the parameters.  Furthermore, in order to decide how to handle
3962
 * the non-integrality, we also need to know whether the coefficients
3963
 * of the other columns in the tableau are integral.  This leads
3964
 * to the following table.  The first two rows do not correspond
3965
 * to a non-integral sample point and are only mentioned for completeness.
3966
 *
3967
 *  constant  parameters  other
3968
 *
3969
 *  int   int   int |
3970
 *  int   int   rat | -> no problem
3971
 *
3972
 *  rat   int   int   -> fail
3973
 *
3974
 *  rat   int   rat   -> cut
3975
 *
3976
 *  int   rat   rat |
3977
 *  rat   rat   rat | -> parametric cut
3978
 *
3979
 *  int   rat   int |
3980
 *  rat   rat   int | -> split context
3981
 *
3982
 * If the parametric constant is completely integral, then there is nothing
3983
 * to be done.  If the constant term is non-integral, but all the other
3984
 * coefficient are integral, then there is nothing that can be done
3985
 * and the tableau has no integral solution.
3986
 * If, on the other hand, one or more of the other columns have rational
3987
 * coefficients, but the parameter coefficients are all integral, then
3988
 * we can perform a regular (non-parametric) cut.
3989
 * Finally, if there is any parameter coefficient that is non-integral,
3990
 * then we need to involve the context tableau.  There are two cases here.
3991
 * If at least one other column has a rational coefficient, then we
3992
 * can perform a parametric cut in the main tableau by adding a new
3993
 * integer division in the context tableau.
3994
 * If all other columns have integral coefficients, then we need to
3995
 * enforce that the rational combination of parameters (c + \sum a_i y_i)/m
3996
 * is always integral.  We do this by introducing an integer division
3997
 * q = floor((c + \sum a_i y_i)/m) and stipulating that its argument should
3998
 * always be integral in the context tableau, i.e., m q = c + \sum a_i y_i.
3999
 * Since q is expressed in the tableau as
4000
 *  c + \sum a_i y_i - m q >= 0
4001
 *  -c - \sum a_i y_i + m q + m - 1 >= 0
4002
 * it is sufficient to add the inequality
4003
 *  -c - \sum a_i y_i + m q >= 0
4004
 * In the part of the context where this inequality does not hold, the
4005
 * main tableau is marked as being empty.
4006
 */
4007
static void find_solutions(struct isl_sol *sol, struct isl_tab *tab)
4008
9.36k
{
4009
9.36k
  struct isl_context *context;
4010
9.36k
  int r;
4011
9.36k
4012
9.36k
  if (
!tab || 9.36k
sol->error9.36k
)
4013
0
    goto error;
4014
9.36k
4015
9.36k
  context = sol->context;
4016
9.36k
4017
9.36k
  if (tab->empty)
4018
0
    goto done;
4019
9.36k
  
if (9.36k
context->op->is_empty(context)9.36k
)
4020
0
    goto done;
4021
9.36k
4022
13.2k
  
for (r = 0; 9.36k
r >= 0 && 13.2k
tab13.2k
&&
!tab->empty13.2k
;
r = restore_lexmin(tab)3.87k
) {
4023
10.6k
    int flags;
4024
10.6k
    int row;
4025
10.6k
    enum isl_tab_row_sign sgn;
4026
10.6k
    int split = -1;
4027
10.6k
    int n_split = 0;
4028
10.6k
4029
72.7k
    for (row = tab->n_redundant; 
row < tab->n_row72.7k
;
++row62.1k
) {
4030
62.1k
      if (!isl_tab_var_from_row(tab, row)->is_nonneg)
4031
705
        continue;
4032
61.4k
      sgn = row_sign(tab, sol, row);
4033
61.4k
      if (!sgn)
4034
0
        goto error;
4035
61.4k
      tab->row_sign[row] = sgn;
4036
61.4k
      if (sgn == isl_tab_row_any)
4037
3.30k
        n_split++;
4038
61.4k
      if (
sgn == isl_tab_row_any && 61.4k
split == -13.30k
)
4039
2.84k
        split = row;
4040
61.4k
      if (sgn == isl_tab_row_neg)
4041
18
        break;
4042
62.1k
    }
4043
10.6k
    
if (10.6k
row < tab->n_row10.6k
)
4044
18
      continue;
4045
10.6k
    
if (10.6k
split != -110.6k
) {
4046
2.84k
      struct isl_vec *ineq;
4047
2.84k
      if (n_split != 1)
4048
331
        split = context->op->best_split(context, tab);
4049
2.84k
      if (split < 0)
4050
0
        goto error;
4051
2.84k
      ineq = get_row_parameter_ineq(tab, split);
4052
2.84k
      if (!ineq)
4053
0
        goto error;
4054
2.84k
      is_strict(ineq);
4055
2.84k
      reset_any_to_unknown(tab);
4056
2.84k
      tab->row_sign[split] = isl_tab_row_pos;
4057
2.84k
      sol_inc_level(sol);
4058
2.84k
      find_in_pos(sol, tab, ineq->el);
4059
2.84k
      tab->row_sign[split] = isl_tab_row_neg;
4060
2.84k
      isl_seq_neg(ineq->el, ineq->el, ineq->size);
4061
2.84k
      isl_int_sub_ui(ineq->el[0], ineq->el[0], 1);
4062
2.84k
      if (!sol->error)
4063
2.84k
        context->op->add_ineq(context, ineq->el, 0, 1);
4064
2.84k
      isl_vec_free(ineq);
4065
2.84k
      if (sol->error)
4066
0
        goto error;
4067
2.84k
      continue;
4068
2.84k
    }
4069
7.77k
    
if (7.77k
tab->rational7.77k
)
4070
3
      break;
4071
7.77k
    row = first_non_integer_row(tab, &flags);
4072
7.77k
    if (row < 0)
4073
6.75k
      break;
4074
1.01k
    
if (1.01k
ISL_FL_ISSET1.01k
(flags, I_PAR)) {
4075
427
      if (
ISL_FL_ISSET427
(flags, I_VAR)) {
4076
0
        if (isl_tab_mark_empty(tab) < 0)
4077
0
          goto error;
4078
0
        break;
4079
0
      }
4080
427
      row = add_cut(tab, row);
4081
1.01k
    } else 
if (592
ISL_FL_ISSET592
(flags, I_VAR)) {
4082
88
      struct isl_vec *div;
4083
88
      struct isl_vec *ineq;
4084
88
      int d;
4085
88
      div = get_row_split_div(tab, row);
4086
88
      if (!div)
4087
0
        goto error;
4088
88
      d = context->op->get_div(context, tab, div);
4089
88
      isl_vec_free(div);
4090
88
      if (d < 0)
4091
0
        goto error;
4092
88
      ineq = ineq_for_div(context->op->peek_basic_set(context), d);
4093
88
      if (!ineq)
4094
0
        goto error;
4095
88
      sol_inc_level(sol);
4096
88
      no_sol_in_strict(sol, tab, ineq);
4097
88
      isl_seq_neg(ineq->el, ineq->el, ineq->size);
4098
88
      context->op->add_ineq(context, ineq->el, 1, 1);
4099
88
      isl_vec_free(ineq);
4100
88
      if (
sol->error || 88
!context->op->is_ok(context)88
)
4101
0
        goto error;
4102
88
      tab = set_row_cst_to_div(tab, row, d);
4103
88
      if (context->op->is_empty(context))
4104
0
        break;
4105
592
    } else
4106
504
      row = add_parametric_cut(tab, row, context);
4107
1.01k
    
if (1.01k
row < 01.01k
)
4108
0
      goto error;
4109
10.6k
  }
4110
9.36k
  
if (9.36k
r < 09.36k
)
4111
0
    goto error;
4112
9.36k
done:
4113
9.36k
  sol_add(sol, tab);
4114
9.36k
  isl_tab_free(tab);
4115
9.36k
  return;
4116
0
error:
4117
0
  isl_tab_free(tab);
4118
0
  sol->error = 1;
4119
0
}
4120
4121
/* Does "sol" contain a pair of partial solutions that could potentially
4122
 * be merged?
4123
 *
4124
 * We currently only check that "sol" is not in an error state
4125
 * and that there are at least two partial solutions of which the final two
4126
 * are defined at the same level.
4127
 */
4128
static int sol_has_mergeable_solutions(struct isl_sol *sol)
4129
6.52k
{
4130
6.52k
  if (sol->error)
4131
0
    return 0;
4132
6.52k
  
if (6.52k
!sol->partial6.52k
)
4133
99
    return 0;
4134
6.42k
  
if (6.42k
!sol->partial->next6.42k
)
4135
4.47k
    return 0;
4136
1.95k
  return sol->partial->level == sol->partial->next->level;
4137
1.95k
}
4138
4139
/* Compute the lexicographic minimum of the set represented by the main
4140
 * tableau "tab" within the context "sol->context_tab".
4141
 *
4142
 * As a preprocessing step, we first transfer all the purely parametric
4143
 * equalities from the main tableau to the context tableau, i.e.,
4144
 * parameters that have been pivoted to a row.
4145
 * These equalities are ignored by the main algorithm, because the
4146
 * corresponding rows may not be marked as being non-negative.
4147
 * In parts of the context where the added equality does not hold,
4148
 * the main tableau is marked as being empty.
4149
 *
4150
 * Before we embark on the actual computation, we save a copy
4151
 * of the context.  When we return, we check if there are any
4152
 * partial solutions that can potentially be merged.  If so,
4153
 * we perform a rollback to the initial state of the context.
4154
 * The merging of partial solutions happens inside calls to
4155
 * sol_dec_level that are pushed onto the undo stack of the context.
4156
 * If there are no partial solutions that can potentially be merged
4157
 * then the rollback is skipped as it would just be wasted effort.
4158
 */
4159
static void find_solutions_main(struct isl_sol *sol, struct isl_tab *tab)
4160
6.52k
{
4161
6.52k
  int row;
4162
6.52k
  void *saved;
4163
6.52k
4164
6.52k
  if (!tab)
4165
0
    goto error;
4166
6.52k
4167
6.52k
  sol->level = 0;
4168
6.52k
4169
72.0k
  for (row = tab->n_redundant; 
row < tab->n_row72.0k
;
++row65.5k
) {
4170
65.5k
    int p;
4171
65.5k
    struct isl_vec *eq;
4172
65.5k
4173
65.5k
    if (tab->row_var[row] < 0)
4174
12.9k
      continue;
4175
52.6k
    
if (52.6k
tab->row_var[row] >= tab->n_param &&
4176
43.5k
        tab->row_var[row] < tab->n_var - tab->n_div)
4177
43.5k
      continue;
4178
9.12k
    
if (9.12k
tab->row_var[row] < tab->n_param9.12k
)
4179
9.12k
      p = tab->row_var[row];
4180
9.12k
    else
4181
0
      p = tab->row_var[row]
4182
0
        + tab->n_param - (tab->n_var - tab->n_div);
4183
9.12k
4184
9.12k
    eq = isl_vec_alloc(tab->mat->ctx, 1+tab->n_param+tab->n_div);
4185
9.12k
    if (!eq)
4186
0
      goto error;
4187
9.12k
    get_row_parameter_line(tab, row, eq->el);
4188
9.12k
    isl_int_neg(eq->el[1 + p], tab->mat->row[row][0]);
4189
9.12k
    eq = isl_vec_normalize(eq);
4190
9.12k
4191
9.12k
    sol_inc_level(sol);
4192
9.12k
    no_sol_in_strict(sol, tab, eq);
4193
9.12k
4194
9.12k
    isl_seq_neg(eq->el, eq->el, eq->size);
4195
9.12k
    sol_inc_level(sol);
4196
9.12k
    no_sol_in_strict(sol, tab, eq);
4197
9.12k
    isl_seq_neg(eq->el, eq->el, eq->size);
4198
9.12k
4199
9.12k
    sol->context->op->add_eq(sol->context, eq->el, 1, 1);
4200
9.12k
4201
9.12k
    isl_vec_free(eq);
4202
9.12k
4203
9.12k
    if (isl_tab_mark_redundant(tab, row) < 0)
4204
0
      goto error;
4205
9.12k
4206
9.12k
    
if (9.12k
sol->context->op->is_empty(sol->context)9.12k
)
4207
0
      break;
4208
9.12k
4209
9.12k
    row = tab->n_redundant - 1;
4210
9.12k
  }
4211
6.52k
4212
6.52k
  saved = sol->context->op->save(sol->context);
4213
6.52k
4214
6.52k
  find_solutions(sol, tab);
4215
6.52k
4216
6.52k
  if (sol_has_mergeable_solutions(sol))
4217
1.87k
    sol->context->op->restore(sol->context, saved);
4218
6.52k
  else
4219
4.65k
    sol->context->op->discard(saved);
4220
6.52k
4221
6.52k
  sol->level = 0;
4222
6.52k
  sol_pop(sol);
4223
6.52k
4224
6.52k
  return;
4225
0
error:
4226
0
  isl_tab_free(tab);
4227
0
  sol->error = 1;
4228
0
}
4229
4230
/* Check if integer division "div" of "dom" also occurs in "bmap".
4231
 * If so, return its position within the divs.
4232
 * If not, return -1.
4233
 */
4234
static int find_context_div(struct isl_basic_map *bmap,
4235
  struct isl_basic_set *dom, unsigned div)
4236
620
{
4237
620
  int i;
4238
620
  unsigned b_dim = isl_space_dim(bmap->dim, isl_dim_all);
4239
620
  unsigned d_dim = isl_space_dim(dom->dim, isl_dim_all);
4240
620
4241
620
  if (isl_int_is_zero(dom->div[div][0]))
4242
130
    return -1;
4243
490
  
if (490
isl_seq_first_non_zero(dom->div[div] + 2 + d_dim, dom->n_div) != -1490
)
4244
2
    return -1;
4245
488
4246
576
  
for (i = 0; 488
i < bmap->n_div576
;
++i88
) {
4247
144
    if (isl_int_is_zero(bmap->div[i][0]))
4248
12
      continue;
4249
132
    
if (132
isl_seq_first_non_zero(bmap->div[i] + 2 + d_dim,
4250
132
             (b_dim - d_dim) + bmap->n_div) != -1)
4251
18
      continue;
4252
114
    
if (114
isl_seq_eq(bmap->div[i], dom->div[div], 2 + d_dim)114
)
4253
56
      return i;
4254
144
  }
4255
432
  return -1;
4256
620
}
4257
4258
/* The correspondence between the variables in the main tableau,
4259
 * the context tableau, and the input map and domain is as follows.
4260
 * The first n_param and the last n_div variables of the main tableau
4261
 * form the variables of the context tableau.
4262
 * In the basic map, these n_param variables correspond to the
4263
 * parameters and the input dimensions.  In the domain, they correspond
4264
 * to the parameters and the set dimensions.
4265
 * The n_div variables correspond to the integer divisions in the domain.
4266
 * To ensure that everything lines up, we may need to copy some of the
4267
 * integer divisions of the domain to the map.  These have to be placed
4268
 * in the same order as those in the context and they have to be placed
4269
 * after any other integer divisions that the map may have.
4270
 * This function performs the required reordering.
4271
 */
4272
static __isl_give isl_basic_map *align_context_divs(
4273
  __isl_take isl_basic_map *bmap, __isl_keep isl_basic_set *dom)
4274
271
{
4275
271
  int i;
4276
271
  int common = 0;
4277
271
  int other;
4278
271
4279
581
  for (i = 0; 
i < dom->n_div581
;
++i310
)
4280
310
    
if (310
find_context_div(bmap, dom, i) != -1310
)
4281
28
      common++;
4282
271
  other = bmap->n_div - common;
4283
271
  if (
dom->n_div - common > 0271
) {
4284
247
    bmap = isl_basic_map_extend_space(bmap, isl_space_copy(bmap->dim),
4285
247
        dom->n_div - common, 0, 0);
4286
247
    if (!bmap)
4287
0
      return NULL;
4288
271
  }
4289
581
  
for (i = 0; 271
i < dom->n_div581
;
++i310
) {
4290
310
    int pos = find_context_div(bmap, dom, i);
4291
310
    if (
pos < 0310
) {
4292
282
      pos = isl_basic_map_alloc_div(bmap);
4293
282
      if (pos < 0)
4294
0
        goto error;
4295
282
      
isl_int_set_si282
(bmap->div[pos][0], 0);
4296
282
    }
4297
310
    
if (310
pos != other + i310
)
4298
15
      isl_basic_map_swap_div(bmap, pos, other + i);
4299
310
  }
4300
271
  return bmap;
4301
0
error:
4302
0
  isl_basic_map_free(bmap);
4303
0
  return NULL;
4304
271
}
4305
4306
/* Base case of isl_tab_basic_map_partial_lexopt, after removing
4307
 * some obvious symmetries.
4308
 *
4309
 * We make sure the divs in the domain are properly ordered,
4310
 * because they will be added one by one in the given order
4311
 * during the construction of the solution map.
4312
 * Furthermore, make sure that the known integer divisions
4313
 * appear before any unknown integer division because the solution
4314
 * may depend on the known integer divisions, while anything that
4315
 * depends on any variable starting from the first unknown integer
4316
 * division is ignored in sol_pma_add.
4317
 */
4318
static struct isl_sol *basic_map_partial_lexopt_base_sol(
4319
  __isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom,
4320
  __isl_give isl_set **empty, int max,
4321
  struct isl_sol *(*init)(__isl_keep isl_basic_map *bmap,
4322
        __isl_take isl_basic_set *dom, int track_empty, int max))
4323
6.85k
{
4324
6.85k
  struct isl_tab *tab;
4325
6.85k
  struct isl_sol *sol = NULL;
4326
6.85k
  struct isl_context *context;
4327
6.85k
4328
6.85k
  if (
dom->n_div6.85k
) {
4329
271
    dom = isl_basic_set_sort_divs(dom);
4330
271
    bmap = align_context_divs(bmap, dom);
4331
271
  }
4332
6.85k
  sol = init(bmap, dom, !!empty, max);
4333
6.85k
  if (!sol)
4334
0
    goto error;
4335
6.85k
4336
6.85k
  context = sol->context;
4337
6.85k
  if (isl_basic_set_plain_is_empty(context->op->peek_basic_set(context)))
4338
12
    /* nothing */;
4339
6.84k
  else 
if (6.84k
isl_basic_map_plain_is_empty(bmap)6.84k
) {
4340
1.03k
    if (sol->add_empty)
4341
1.03k
      sol->add_empty(sol,
4342
1.03k
        isl_basic_set_copy(context->op->peek_basic_set(context)));
4343
6.84k
  } else {
4344
5.80k
    tab = tab_for_lexmin(bmap,
4345
5.80k
            context->op->peek_basic_set(context), 1, max);
4346
5.80k
    tab = context->op->detect_nonnegative_parameters(context, tab);
4347
5.80k
    find_solutions_main(sol, tab);
4348
5.80k
  }
4349
6.85k
  if (sol->error)
4350
0
    goto error;
4351
6.85k
4352
6.85k
  isl_basic_map_free(bmap);
4353
6.85k
  return sol;
4354
0
error:
4355
0
  sol_free(sol);
4356
0
  isl_basic_map_free(bmap);
4357
0
  return NULL;
4358
6.85k
}
4359
4360
/* Base case of isl_tab_basic_map_partial_lexopt, after removing
4361
 * some obvious symmetries.
4362
 *
4363
 * We call basic_map_partial_lexopt_base_sol and extract the results.
4364
 */
4365
static __isl_give isl_map *basic_map_partial_lexopt_base(
4366
  __isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom,
4367
  __isl_give isl_set **empty, int max)
4368
3.40k
{
4369
3.40k
  isl_map *result = NULL;
4370
3.40k
  struct isl_sol *sol;
4371
3.40k
  struct isl_sol_map *sol_map;
4372
3.40k
4373
3.40k
  sol = basic_map_partial_lexopt_base_sol(bmap, dom, empty, max,
4374
3.40k
            &sol_map_init);
4375
3.40k
  if (!sol)
4376
0
    return NULL;
4377
3.40k
  sol_map = (struct isl_sol_map *) sol;
4378
3.40k
4379
3.40k
  result = isl_map_copy(sol_map->map);
4380
3.40k
  if (empty)
4381
2.60k
    *empty = isl_set_copy(sol_map->empty);
4382
3.40k
  sol_free(&sol_map->sol);
4383
3.40k
  return result;
4384
3.40k
}
4385
4386
/* Return a count of the number of occurrences of the "n" first
4387
 * variables in the inequality constraints of "bmap".
4388
 */
4389
static __isl_give int *count_occurrences(__isl_keep isl_basic_map *bmap,
4390
  int n)
4391
7.01k
{
4392
7.01k
  int i, j;
4393
7.01k
  isl_ctx *ctx;
4394
7.01k
  int *occurrences;
4395
7.01k
4396
7.01k
  if (!bmap)
4397
0
    return NULL;
4398
7.01k
  ctx = isl_basic_map_get_ctx(bmap);
4399
7.01k
  occurrences = isl_calloc_array(ctx, int, n);
4400
7.01k
  if (!occurrences)
4401
0
    return NULL;
4402
7.01k
4403
23.9k
  
for (i = 0; 7.01k
i < bmap->n_ineq23.9k
;
++i16.9k
) {
4404
90.4k
    for (j = 0; 
j < n90.4k
;
++j73.4k
) {
4405
73.4k
      if (
!73.4k
isl_int_is_zero73.4k
(bmap->ineq[i][1 + j]))
4406
13.6k
        occurrences[j]++;
4407
73.4k
    }
4408
16.9k
  }
4409
7.01k
4410
7.01k
  return occurrences;
4411
7.01k
}
4412
4413
/* Do all of the "n" variables with non-zero coefficients in "c"
4414
 * occur in exactly a single constraint.
4415
 * "occurrences" is an array of length "n" containing the number
4416
 * of occurrences of each of the variables in the inequality constraints.
4417
 */
4418
static int single_occurrence(int n, isl_int *c, int *occurrences)
4419
7.84k
{
4420
7.84k
  int i;
4421
7.84k
4422
28.6k
  for (i = 0; 
i < n28.6k
;
++i20.8k
) {
4423
22.9k
    if (isl_int_is_zero(c[i]))
4424
20.1k
      continue;
4425
2.84k
    
if (2.84k
occurrences[i] != 12.84k
)
4426
2.15k
      return 0;
4427
22.9k
  }
4428
7.84k
4429
5.69k
  return 1;
4430
7.84k
}
4431
4432
/* Do all of the "n" initial variables that occur in inequality constraint
4433
 * "ineq" of "bmap" only occur in that constraint?
4434
 */
4435
static int all_single_occurrence(__isl_keep isl_basic_map *bmap, int ineq,
4436
  int n)
4437
3
{
4438
3
  int i, j;
4439
3
4440
21
  for (i = 0; 
i < n21
;
++i18
) {
4441
18
    if (isl_int_is_zero(bmap->ineq[ineq][1 + i]))
4442
16
      continue;
4443
10
    
for (j = 0; 2
j < bmap->n_ineq10
;
++j8
) {
4444
8
      if (j == ineq)
4445
2
        continue;
4446
6
      
if (6
!6
isl_int_is_zero6
(bmap->ineq[j][1 + i]))
4447
0
        return 0;
4448
8
    }
4449
18
  }
4450
3
4451
3
  return 1;
4452
3
}
4453
4454
/* Structure used during detection of parallel constraints.
4455
 * n_in: number of "input" variables: isl_dim_param + isl_dim_in
4456
 * n_out: number of "output" variables: isl_dim_out + isl_dim_div
4457
 * val: the coefficients of the output variables
4458
 */
4459
struct isl_constraint_equal_info {
4460
  unsigned n_in;
4461
  unsigned n_out;
4462
  isl_int *val;
4463
};
4464
4465
/* Check whether the coefficients of the output variables
4466
 * of the constraint in "entry" are equal to info->val.
4467
 */
4468
static int constraint_equal(const void *entry, const void *val)
4469
161
{
4470
161
  isl_int **row = (isl_int **)entry;
4471
161
  const struct isl_constraint_equal_info *info = val;
4472
161
4473
161
  return isl_seq_eq((*row) + 1 + info->n_in, info->val, info->n_out);
4474
161
}
4475
4476
/* Check whether "bmap" has a pair of constraints that have
4477
 * the same coefficients for the output variables.
4478
 * Note that the coefficients of the existentially quantified
4479
 * variables need to be zero since the existentially quantified
4480
 * of the result are usually not the same as those of the input.
4481
 * Furthermore, check that each of the input variables that occur
4482
 * in those constraints does not occur in any other constraint.
4483
 * If so, return true and return the row indices of the two constraints
4484
 * in *first and *second.
4485
 */
4486
static isl_bool parallel_constraints(__isl_keep isl_basic_map *bmap,
4487
  int *first, int *second)
4488
7.01k
{
4489
7.01k
  int i;
4490
7.01k
  isl_ctx *ctx;
4491
7.01k
  int *occurrences = NULL;
4492
7.01k
  struct isl_hash_table *table = NULL;
4493
7.01k
  struct isl_hash_table_entry *entry;
4494
7.01k
  struct isl_constraint_equal_info info;
4495
7.01k
  unsigned n_out;
4496
7.01k
  unsigned n_div;
4497
7.01k
4498
7.01k
  ctx = isl_basic_map_get_ctx(bmap);
4499
7.01k
  table = isl_hash_table_alloc(ctx, bmap->n_ineq);
4500
7.01k
  if (!table)
4501
0
    goto error;
4502
7.01k
4503
7.01k
  info.n_in = isl_basic_map_dim(bmap, isl_dim_param) +
4504
7.01k
        isl_basic_map_dim(bmap, isl_dim_in);
4505
7.01k
  occurrences = count_occurrences(bmap, info.n_in);
4506
7.01k
  if (
info.n_in && 7.01k
!occurrences6.18k
)
4507
0
    goto error;
4508
7.01k
  n_out = isl_basic_map_dim(bmap, isl_dim_out);
4509
7.01k
  n_div = isl_basic_map_dim(bmap, isl_dim_div);
4510
7.01k
  info.n_out = n_out + n_div;
4511
23.6k
  for (i = 0; 
i < bmap->n_ineq23.6k
;
++i16.6k
) {
4512
16.8k
    uint32_t hash;
4513
16.8k
4514
16.8k
    info.val = bmap->ineq[i] + 1 + info.n_in;
4515
16.8k
    if (isl_seq_first_non_zero(info.val, n_out) < 0)
4516
8.88k
      continue;
4517
7.94k
    
if (7.94k
isl_seq_first_non_zero(info.val + n_out, n_div) >= 07.94k
)
4518
94
      continue;
4519
7.84k
    
if (7.84k
!single_occurrence(info.n_in, bmap->ineq[i] + 1,
4520
7.84k
          occurrences))
4521
2.15k
      continue;
4522
5.69k
    hash = isl_seq_get_hash(info.val, info.n_out);
4523
5.69k
    entry = isl_hash_table_find(ctx, table, hash,
4524
5.69k
              constraint_equal, &info, 1);
4525
5.69k
    if (!entry)
4526
0
      goto error;
4527
5.69k
    
if (5.69k
entry->data5.69k
)
4528
161
      break;
4529
5.53k
    entry->data = &bmap->ineq[i];
4530
5.53k
  }
4531
7.01k
4532
7.01k
  
if (7.01k
i < bmap->n_ineq7.01k
) {
4533
161
    *first = ((isl_int **)entry->data) - bmap->ineq; 
4534
161
    *second = i;
4535
161
  }
4536
7.01k
4537
7.01k
  isl_hash_table_free(ctx, table);
4538
7.01k
  free(occurrences);
4539
7.01k
4540
7.01k
  return i < bmap->n_ineq;
4541
0
error:
4542
0
  isl_hash_table_free(ctx, table);
4543
0
  free(occurrences);
4544
0
  return isl_bool_error;
4545
7.01k
}
4546
4547
/* Given a set of upper bounds in "var", add constraints to "bset"
4548
 * that make the i-th bound smallest.
4549
 *
4550
 * In particular, if there are n bounds b_i, then add the constraints
4551
 *
4552
 *  b_i <= b_j  for j > i
4553
 *  b_i <  b_j  for j < i
4554
 */
4555
static __isl_give isl_basic_set *select_minimum(__isl_take isl_basic_set *bset,
4556
  __isl_keep isl_mat *var, int i)
4557
548
{
4558
548
  isl_ctx *ctx;
4559
548
  int j, k;
4560
548
4561
548
  ctx = isl_mat_get_ctx(var);
4562
548
4563
1.66k
  for (j = 0; 
j < var->n_row1.66k
;
++j1.11k
) {
4564
1.11k
    if (j == i)
4565
548
      continue;
4566
566
    k = isl_basic_set_alloc_inequality(bset);
4567
566
    if (k < 0)
4568
0
      goto error;
4569
566
    isl_seq_combine(bset->ineq[k], ctx->one, var->row[j],
4570
566
        ctx->negone, var->row[i], var->n_col);
4571
566
    isl_int_set_si(bset->ineq[k][var->n_col], 0);
4572
566
    if (j < i)
4573
283
      isl_int_sub_ui(bset->ineq[k][0], bset->ineq[k][0], 1);
4574
1.11k
  }
4575
548
4576
548
  bset = isl_basic_set_finalize(bset);
4577
548
4578
548
  return bset;
4579
0
error:
4580
0
  isl_basic_set_free(bset);
4581
0
  return NULL;
4582
548
}
4583
4584
/* Given a set of upper bounds on the last "input" variable m,
4585
 * construct a set that assigns the minimal upper bound to m, i.e.,
4586
 * construct a set that divides the space into cells where one
4587
 * of the upper bounds is smaller than all the others and assign
4588
 * this upper bound to m.
4589
 *
4590
 * In particular, if there are n bounds b_i, then the result
4591
 * consists of n basic sets, each one of the form
4592
 *
4593
 *  m = b_i
4594
 *  b_i <= b_j  for j > i
4595
 *  b_i <  b_j  for j < i
4596
 */
4597
static __isl_give isl_set *set_minimum(__isl_take isl_space *dim,
4598
  __isl_take isl_mat *var)
4599
161
{
4600
161
  int i, k;
4601
161
  isl_basic_set *bset = NULL;
4602
161
  isl_set *set = NULL;
4603
161
4604
161
  if (
!dim || 161
!var161
)
4605
0
    goto error;
4606
161
4607
161
  set = isl_set_alloc_space(isl_space_copy(dim),
4608
161
        var->n_row, ISL_SET_DISJOINT);
4609
161
4610
486
  for (i = 0; 
i < var->n_row486
;
++i325
) {
4611
325
    bset = isl_basic_set_alloc_space(isl_space_copy(dim), 0,
4612
325
                 1, var->n_row - 1);
4613
325
    k = isl_basic_set_alloc_equality(bset);
4614
325
    if (k < 0)
4615
0
      goto error;
4616
325
    isl_seq_cpy(bset->eq[k], var->row[i], var->n_col);
4617
325
    isl_int_set_si(bset->eq[k][var->n_col], -1);
4618
325
    bset = select_minimum(bset, var, i);
4619
325
    set = isl_set_add_basic_set(set, bset);
4620
325
  }
4621
161
4622
161
  isl_space_free(dim);
4623
161
  isl_mat_free(var);
4624
161
  return set;
4625
0
error:
4626
0
  isl_basic_set_free(bset);
4627
0
  isl_set_free(set);
4628
0
  isl_space_free(dim);
4629
0
  isl_mat_free(var);
4630
0
  return NULL;
4631
161
}
4632
4633
/* Given that the last input variable of "bmap" represents the minimum
4634
 * of the bounds in "cst", check whether we need to split the domain
4635
 * based on which bound attains the minimum.
4636
 *
4637
 * A split is needed when the minimum appears in an integer division
4638
 * or in an equality.  Otherwise, it is only needed if it appears in
4639
 * an upper bound that is different from the upper bounds on which it
4640
 * is defined.
4641
 */
4642
static isl_bool need_split_basic_map(__isl_keep isl_basic_map *bmap,
4643
  __isl_keep isl_mat *cst)
4644
111
{
4645
111
  int i, j;
4646
111
  unsigned total;
4647
111
  unsigned pos;
4648
111
4649
111
  pos = cst->n_col - 1;
4650
111
  total = isl_basic_map_dim(bmap, isl_dim_all);
4651
111
4652
145
  for (i = 0; 
i < bmap->n_div145
;
++i34
)
4653
108
    
if (108
!108
isl_int_is_zero108
(bmap->div[i][2 + pos]))
4654
74
      return isl_bool_true;
4655
111
4656
79
  
for (i = 0; 37
i < bmap->n_eq79
;
++i42
)
4657
46
    
if (46
!46
isl_int_is_zero46
(bmap->eq[i][1 + pos]))
4658
4
      return isl_bool_true;
4659
37
4660
149
  
for (i = 0; 33
i < bmap->n_ineq149
;
++i116
) {
4661
136
    if (isl_int_is_nonneg(bmap->ineq[i][1 + pos]))
4662
78
      continue;
4663
58
    
if (58
!58
isl_int_is_negone58
(bmap->ineq[i][1 + pos]))
4664
0
      return isl_bool_true;
4665
58
    
if (58
isl_seq_first_non_zero(bmap->ineq[i] + 1 + pos + 1,
4666
58
             total - pos - 1) >= 0)
4667
14
      return isl_bool_true;
4668
44
4669
75
    
for (j = 0; 44
j < cst->n_row75
;
++j31
)
4670
69
      
if (69
isl_seq_eq(bmap->ineq[i], cst->row[j], cst->n_col)69
)
4671
38
        break;
4672
44
    if (j >= cst->n_row)
4673
6
      return isl_bool_true;
4674
136
  }
4675
33
4676
13
  return isl_bool_false;
4677
111
}
4678
4679
/* Given that the last set variable of "bset" represents the minimum
4680
 * of the bounds in "cst", check whether we need to split the domain
4681
 * based on which bound attains the minimum.
4682
 *
4683
 * We simply call need_split_basic_map here.  This is safe because
4684
 * the position of the minimum is computed from "cst" and not
4685
 * from "bmap".
4686
 */
4687
static isl_bool need_split_basic_set(__isl_keep isl_basic_set *bset,
4688
  __isl_keep isl_mat *cst)
4689
13
{
4690
13
  return need_split_basic_map(bset_to_bmap(bset), cst);
4691
13
}
4692
4693
/* Given that the last set variable of "set" represents the minimum
4694
 * of the bounds in "cst", check whether we need to split the domain
4695
 * based on which bound attains the minimum.
4696
 */
4697
static isl_bool need_split_set(__isl_keep isl_set *set, __isl_keep isl_mat *cst)
4698
7
{
4699
7
  int i;
4700
7
4701
14
  for (i = 0; 
i < set->n14
;
++i7
) {
4702
7
    isl_bool split;
4703
7
4704
7
    split = need_split_basic_set(set->p[i], cst);
4705
7
    if (
split < 0 || 7
split7
)
4706
0
      return split;
4707
7
  }
4708
7
4709
7
  return isl_bool_false;
4710
7
}
4711
4712
/* Given a set of which the last set variable is the minimum
4713
 * of the bounds in "cst", split each basic set in the set
4714
 * in pieces where one of the bounds is (strictly) smaller than the others.
4715
 * This subdivision is given in "min_expr".
4716
 * The variable is subsequently projected out.
4717
 *
4718
 * We only do the split when it is needed.
4719
 * For example if the last input variable m = min(a,b) and the only
4720
 * constraints in the given basic set are lower bounds on m,
4721
 * i.e., l <= m = min(a,b), then we can simply project out m
4722
 * to obtain l <= a and l <= b, without having to split on whether
4723
 * m is equal to a or b.
4724
 */
4725
static __isl_give isl_set *split(__isl_take isl_set *empty,
4726
  __isl_take isl_set *min_expr, __isl_take isl_mat *cst)
4727
6
{
4728
6
  int n_in;
4729
6
  int i;
4730
6
  isl_space *dim;
4731
6
  isl_set *res;
4732
6
4733
6
  if (
!empty || 6
!min_expr6
||
!cst6
)
4734
0
    goto error;
4735
6
4736
6
  n_in = isl_set_dim(empty, isl_dim_set);
4737
6
  dim = isl_set_get_space(empty);
4738
6
  dim = isl_space_drop_dims(dim, isl_dim_set, n_in - 1, 1);
4739
6
  res = isl_set_empty(dim);
4740
6
4741
12
  for (i = 0; 
i < empty->n12
;
++i6
) {
4742
6
    isl_bool split;
4743
6
    isl_set *set;
4744
6
4745
6
    set = isl_set_from_basic_set(isl_basic_set_copy(empty->p[i]));
4746
6
    split = need_split_basic_set(empty->p[i], cst);
4747
6
    if (split < 0)
4748
0
      set = isl_set_free(set);
4749
6
    else 
if (6
split6
)
4750
6
      set = isl_set_intersect(set, isl_set_copy(min_expr));
4751
6
    set = isl_set_remove_dims(set, isl_dim_set, n_in - 1, 1);
4752
6
4753
6
    res = isl_set_union_disjoint(res, set);
4754
6
  }
4755
6
4756
6
  isl_set_free(empty);
4757
6
  isl_set_free(min_expr);
4758
6
  isl_mat_free(cst);
4759
6
  return res;
4760
0
error:
4761
0
  isl_set_free(empty);
4762
0
  isl_set_free(min_expr);
4763
0
  isl_mat_free(cst);
4764
0
  return NULL;
4765
6
}
4766
4767
/* Given a map of which the last input variable is the minimum
4768
 * of the bounds in "cst", split each basic set in the set
4769
 * in pieces where one of the bounds is (strictly) smaller than the others.
4770
 * This subdivision is given in "min_expr".
4771
 * The variable is subsequently projected out.
4772
 *
4773
 * The implementation is essentially the same as that of "split".
4774
 */
4775
static __isl_give isl_map *split_domain(__isl_take isl_map *opt,
4776
  __isl_take isl_set *min_expr, __isl_take isl_mat *cst)
4777
51
{
4778
51
  int n_in;
4779
51
  int i;
4780
51
  isl_space *dim;
4781
51
  isl_map *res;
4782
51
4783
51
  if (
!opt || 51
!min_expr51
||
!cst51
)
4784
0
    goto error;
4785
51
4786
51
  n_in = isl_map_dim(opt, isl_dim_in);
4787
51
  dim = isl_map_get_space(opt);
4788
51
  dim = isl_space_drop_dims(dim, isl_dim_in, n_in - 1, 1);
4789
51
  res = isl_map_empty(dim);
4790
51
4791
149
  for (i = 0; 
i < opt->n149
;
++i98
) {
4792
98
    isl_map *map;
4793
98
    isl_bool split;
4794
98
4795
98
    map = isl_map_from_basic_map(isl_basic_map_copy(opt->p[i]));
4796
98
    split = need_split_basic_map(opt->p[i], cst);
4797
98
    if (split < 0)
4798
0
      map = isl_map_free(map);
4799
98
    else 
if (98
split98
)
4800
92
      map = isl_map_intersect_domain(map,
4801
92
                   isl_set_copy(min_expr));
4802
98
    map = isl_map_remove_dims(map, isl_dim_in, n_in - 1, 1);
4803
98
4804
98
    res = isl_map_union_disjoint(res, map);
4805
98
  }
4806
51
4807
51
  isl_map_free(opt);
4808
51
  isl_set_free(min_expr);
4809
51
  isl_mat_free(cst);
4810
51
  return res;
4811
0
error:
4812
0
  isl_map_free(opt);
4813
0
  isl_set_free(min_expr);
4814
0
  isl_mat_free(cst);
4815
0
  return NULL;
4816
51
}
4817
4818
static __isl_give isl_map *basic_map_partial_lexopt(
4819
  __isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom,
4820
  __isl_give isl_set **empty, int max);
4821
4822
/* This function is called from basic_map_partial_lexopt_symm.
4823
 * The last variable of "bmap" and "dom" corresponds to the minimum
4824
 * of the bounds in "cst".  "map_space" is the space of the original
4825
 * input relation (of basic_map_partial_lexopt_symm) and "set_space"
4826
 * is the space of the original domain.
4827
 *
4828
 * We recursively call basic_map_partial_lexopt and then plug in
4829
 * the definition of the minimum in the result.
4830
 */
4831
static __isl_give isl_map *basic_map_partial_lexopt_symm_core(
4832
  __isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom,
4833
  __isl_give isl_set **empty, int max, __isl_take isl_mat *cst,
4834
  __isl_take isl_space *map_space, __isl_take isl_space *set_space)
4835
51
{
4836
51
  isl_map *opt;
4837
51
  isl_set *min_expr;
4838
51
4839
51
  min_expr = set_minimum(isl_basic_set_get_space(dom), isl_mat_copy(cst));
4840
51
4841
51
  opt = basic_map_partial_lexopt(bmap, dom, empty, max);
4842
51
4843
51
  if (
empty51
) {
4844
6
    *empty = split(*empty,
4845
6
             isl_set_copy(min_expr), isl_mat_copy(cst));
4846
6
    *empty = isl_set_reset_space(*empty, set_space);
4847
6
  }
4848
51
4849
51
  opt = split_domain(opt, min_expr, cst);
4850
51
  opt = isl_map_reset_space(opt, map_space);
4851
51
4852
51
  return opt;
4853
51
}
4854
4855
/* Extract a domain from "bmap" for the purpose of computing
4856
 * a lexicographic optimum.
4857
 *
4858
 * This function is only called when the caller wants to compute a full
4859
 * lexicographic optimum, i.e., without specifying a domain.  In this case,
4860
 * the caller is not interested in the part of the domain space where
4861
 * there is no solution and the domain can be initialized to those constraints
4862
 * of "bmap" that only involve the parameters and the input dimensions.
4863
 * This relieves the parametric programming engine from detecting those
4864
 * inequalities and transferring them to the context.  More importantly,
4865
 * it ensures that those inequalities are transferred first and not
4866
 * intermixed with inequalities that actually split the domain.
4867
 *
4868
 * If the caller does not require the absence of existentially quantified
4869
 * variables in the result (i.e., if ISL_OPT_QE is not set in "flags"),
4870
 * then the actual domain of "bmap" can be used.  This ensures that
4871
 * the domain does not need to be split at all just to separate out
4872
 * pieces of the domain that do not have a solution from piece that do.
4873
 * This domain cannot be used in general because it may involve
4874
 * (unknown) existentially quantified variables which will then also
4875
 * appear in the solution.
4876
 */
4877
static __isl_give isl_basic_set *extract_domain(__isl_keep isl_basic_map *bmap,
4878
  unsigned flags)
4879
3.16k
{
4880
3.16k
  int n_div;
4881
3.16k
  int n_out;
4882
3.16k
4883
3.16k
  n_div = isl_basic_map_dim(bmap, isl_dim_div);
4884
3.16k
  n_out = isl_basic_map_dim(bmap, isl_dim_out);
4885
3.16k
  bmap = isl_basic_map_copy(bmap);
4886
3.16k
  if (
ISL_FL_ISSET3.16k
(flags, ISL_OPT_QE)) {
4887
168
    bmap = isl_basic_map_drop_constraints_involving_dims(bmap,
4888
168
              isl_dim_div, 0, n_div);
4889
168
    bmap = isl_basic_map_drop_constraints_involving_dims(bmap,
4890
168
              isl_dim_out, 0, n_out);
4891
168
  }
4892
3.16k
  return isl_basic_map_domain(bmap);
4893
3.16k
}
4894
4895
#undef TYPE
4896
#define TYPE  isl_map
4897
#undef SUFFIX
4898
#define SUFFIX
4899
#include "isl_tab_lexopt_templ.c"
4900
4901
struct isl_sol_for {
4902
  struct isl_sol  sol;
4903
  isl_stat  (*fn)(__isl_take isl_basic_set *dom,
4904
        __isl_take isl_aff_list *list, void *user);
4905
  void    *user;
4906
};
4907
4908
static void sol_for_free(struct isl_sol *sol)
4909
722
{
4910
722
}
4911
4912
/* Add the solution identified by the tableau and the context tableau.
4913
 * In particular, "dom" represents the context and "ma" expresses
4914
 * the solution on that context.
4915
 *
4916
 * See documentation of sol_add for more details.
4917
 *
4918
 * Instead of constructing a basic map, this function calls a user
4919
 * defined function with the current context as a basic set and
4920
 * a list of affine expressions representing the relation between
4921
 * the input and output.  The space over which the affine expressions
4922
 * are defined is the same as that of the domain.  The number of
4923
 * affine expressions in the list is equal to the number of output variables.
4924
 */
4925
static void sol_for_add(struct isl_sol_for *sol,
4926
  __isl_take isl_basic_set *dom, __isl_take isl_multi_aff *ma)
4927
723
{
4928
723
  int i, n;
4929
723
  isl_ctx *ctx;
4930
723
  isl_aff *aff;
4931
723
  isl_aff_list *list;
4932
723
4933
723
  if (
sol->sol.error || 723
!dom723
||
!ma723
)
4934
0
    goto error;
4935
723
4936
723
  ctx = isl_basic_set_get_ctx(dom);
4937
723
  n = isl_multi_aff_dim(ma, isl_dim_out);
4938
723
  list = isl_aff_list_alloc(ctx, n);
4939
1.44k
  for (i = 0; 
i < n1.44k
;
++i723
) {
4940
723
    aff = isl_multi_aff_get_aff(ma, i);
4941
723
    list = isl_aff_list_add(list, aff);
4942
723
  }
4943
723
4944
723
  dom = isl_basic_set_finalize(dom);
4945
723
4946
723
  if (sol->fn(isl_basic_set_copy(dom), list, sol->user) < 0)
4947
0
    goto error;
4948
723
4949
723
  isl_basic_set_free(dom);
4950
723
  isl_multi_aff_free(ma);
4951
723
  return;
4952
0
error:
4953
0
  isl_basic_set_free(dom);
4954
0
  isl_multi_aff_free(ma);
4955
0
  sol->sol.error = 1;
4956
0
}
4957
4958
static void sol_for_add_wrap(struct isl_sol *sol,
4959
  __isl_take isl_basic_set *dom, __isl_take isl_multi_aff *ma)
4960
723
{
4961
723
  sol_for_add((struct isl_sol_for *)sol, dom, ma);
4962
723
}
4963
4964
static struct isl_sol_for *sol_for_init(__isl_keep isl_basic_map *bmap, int max,
4965
  isl_stat (*fn)(__isl_take isl_basic_set *dom,
4966
    __isl_take isl_aff_list *list, void *user),
4967
  void *user)
4968
722
{
4969
722
  struct isl_sol_for *sol_for = NULL;
4970
722
  isl_space *dom_dim;
4971
722
  struct isl_basic_set *dom = NULL;
4972
722
4973
722
  sol_for = isl_calloc_type(bmap->ctx, struct isl_sol_for);
4974
722
  if (!sol_for)
4975
0
    goto error;
4976
722
4977
722
  dom_dim = isl_space_domain(isl_space_copy(bmap->dim));
4978
722
  dom = isl_basic_set_universe(dom_dim);
4979
722
4980
722
  sol_for->sol.free = &sol_for_free;
4981
722
  if (sol_init(&sol_for->sol, bmap, dom, max) < 0)
4982
0
    goto error;
4983
722
  sol_for->fn = fn;
4984
722
  sol_for->user = user;
4985
722
  sol_for->sol.add = &sol_for_add_wrap;
4986
722
  sol_for->sol.add_empty = NULL;
4987
722
4988
722
  isl_basic_set_free(dom);
4989
722
  return sol_for;
4990
0
error:
4991
0
  isl_basic_set_free(dom);
4992
0
  sol_free(&sol_for->sol);
4993
0
  return NULL;
4994
722
}
4995
4996
static void sol_for_find_solutions(struct isl_sol_for *sol_for,
4997
  struct isl_tab *tab)
4998
722
{
4999
722
  find_solutions_main(&sol_for->sol, tab);
5000
722
}
5001
5002
isl_stat isl_basic_map_foreach_lexopt(__isl_keep isl_basic_map *bmap, int max,
5003
  isl_stat (*fn)(__isl_take isl_basic_set *dom,
5004
    __isl_take isl_aff_list *list, void *user),
5005
  void *user)
5006
722
{
5007
722
  struct isl_sol_for *sol_for = NULL;
5008
722
5009
722
  bmap = isl_basic_map_copy(bmap);
5010
722
  bmap = isl_basic_map_detect_equalities(bmap);
5011
722
  if (!bmap)
5012
0
    return isl_stat_error;
5013
722
5014
722
  sol_for = sol_for_init(bmap, max, fn, user);
5015
722
  if (!sol_for)
5016
0
    goto error;
5017
722
5018
722
  
if (722
isl_basic_map_plain_is_empty(bmap)722
)
5019
0
    /* nothing */;
5020
722
  else {
5021
722
    struct isl_tab *tab;
5022
722
    struct isl_context *context = sol_for->sol.context;
5023
722
    tab = tab_for_lexmin(bmap,
5024
722
        context->op->peek_basic_set(context), 1, max);
5025
722
    tab = context->op->detect_nonnegative_parameters(context, tab);
5026
722
    sol_for_find_solutions(sol_for, tab);
5027
722
    if (sol_for->sol.error)
5028
0
      goto error;
5029
722
  }
5030
722
5031
722
  sol_free(&sol_for->sol);
5032
722
  isl_basic_map_free(bmap);
5033
722
  return isl_stat_ok;
5034
0
error:
5035
0
  sol_free(&sol_for->sol);
5036
0
  isl_basic_map_free(bmap);
5037
0
  return isl_stat_error;
5038
722
}
5039
5040
/* Extract the subsequence of the sample value of "tab"
5041
 * starting at "pos" and of length "len".
5042
 */
5043
static __isl_give isl_vec *extract_sample_sequence(struct isl_tab *tab,
5044
  int pos, int len)
5045
797
{
5046
797
  int i;
5047
797
  isl_ctx *ctx;
5048
797
  isl_vec *v;
5049
797
5050
797
  ctx = isl_tab_get_ctx(tab);
5051
797
  v = isl_vec_alloc(ctx, len);
5052
797
  if (!v)
5053
0
    return NULL;
5054
4.11k
  
for (i = 0; 797
i < len4.11k
;
++i3.32k
) {
5055
3.32k
    if (
!tab->var[pos + i].is_row3.32k
) {
5056
2.52k
      isl_int_set_si(v->el[i], 0);
5057
3.32k
    } else {
5058
793
      int row;
5059
793
5060
793
      row = tab->var[pos + i].index;
5061
793
      isl_int_set(v->el[i], tab->mat->row[row][1]);
5062
793
    }
5063
3.32k
  }
5064
797
5065
797
  return v;
5066
797
}
5067
5068
/* Check if the sequence of variables starting at "pos"
5069
 * represents a trivial solution according to "trivial".
5070
 * That is, is the result of applying "trivial" to this sequence
5071
 * equal to the zero vector?
5072
 */
5073
static isl_bool region_is_trivial(struct isl_tab *tab, int pos,
5074
  __isl_keep isl_mat *trivial)
5075
909
{
5076
909
  int n, len;
5077
909
  isl_vec *v;
5078
909
  isl_bool is_trivial;
5079
909
5080
909
  if (!trivial)
5081
0
    return isl_bool_error;
5082
909
5083
909
  n = isl_mat_rows(trivial);
5084
909
  if (n == 0)
5085
112
    return isl_bool_false;
5086
797
5087
797
  len = isl_mat_cols(trivial);
5088
797
  v = extract_sample_sequence(tab, pos, len);
5089
797
  v = isl_mat_vec_product(isl_mat_copy(trivial), v);
5090
797
  is_trivial = isl_vec_is_zero(v);
5091
797
  isl_vec_free(v);
5092
797
5093
797
  return is_trivial;
5094
797
}
5095
5096
/* Return the index of the first trivial region, "n_region" if all regions
5097
 * are non-trivial or -1 in case of error.
5098
 */
5099
static int first_trivial_region(struct isl_tab *tab,
5100
  int n_region, struct isl_trivial_region *region)
5101
721
{
5102
721
  int i;
5103
721
5104
1.20k
  for (i = 0; 
i < n_region1.20k
;
++i486
) {
5105
909
    isl_bool trivial;
5106
909
    trivial = region_is_trivial(tab, region[i].pos,
5107
909
          region[i].trivial);
5108
909
    if (trivial < 0)
5109
0
      return -1;
5110
909
    
if (909
trivial909
)
5111
423
      return i;
5112
909
  }
5113
721
5114
298
  return n_region;
5115
721
}
5116
5117
/* Check if the solution is optimal, i.e., whether the first
5118
 * n_op entries are zero.
5119
 */
5120
static int is_optimal(__isl_keep isl_vec *sol, int n_op)
5121
298
{
5122
298
  int i;
5123
298
5124
776
  for (i = 0; 
i < n_op776
;
++i478
)
5125
590
    
if (590
!590
isl_int_is_zero590
(sol->el[1 + i]))
5126
112
      return 0;
5127
186
  return 1;
5128
298
}
5129
5130
/* Add constraints to "tab" that ensure that any solution is significantly
5131
 * better than that represented by "sol".  That is, find the first
5132
 * relevant (within first n_op) non-zero coefficient and force it (along
5133
 * with all previous coefficients) to be zero.
5134
 * If the solution is already optimal (all relevant coefficients are zero),
5135
 * then just mark the table as empty.
5136
 * "n_zero" is the number of coefficients that have been forced zero
5137
 * by previous calls to this function at the same level.
5138
 * Return the updated number of forced zero coefficients or -1 on error.
5139
 *
5140
 * This function assumes that at least 2 * (n_op - n_zero) more rows and
5141
 * at least 2 * (n_op - n_zero) more elements in the constraint array
5142
 * are available in the tableau.
5143
 */
5144
static int force_better_solution(struct isl_tab *tab,
5145
  __isl_keep isl_vec *sol, int n_op, int n_zero)
5146
123
{
5147
123
  int i, n;
5148
123
  isl_ctx *ctx;
5149
123
  isl_vec *v = NULL;
5150
123
5151
123
  if (!sol)
5152
0
    return -1;
5153
123
5154
237
  
for (i = n_zero; 123
i < n_op237
;
++i114
)
5155
237
    
if (237
!237
isl_int_is_zero237
(sol->el[1 + i]))
5156
123
      break;
5157
123
5158
123
  if (
i == n_op123
) {
5159
0
    if (isl_tab_mark_empty(tab) < 0)
5160
0
      return -1;
5161
0
    return n_op;
5162
0
  }
5163
123
5164
123
  ctx = isl_vec_get_ctx(sol);
5165
123
  v = isl_vec_alloc(ctx, 1 + tab->n_var);
5166
123
  if (!v)
5167
0
    return -1;
5168
123
5169
123
  n = i + 1;
5170
360
  for (; 
i >= n_zero360
;
--i237
) {
5171
237
    v = isl_vec_clr(v);
5172
237
    isl_int_set_si(v->el[1 + i], -1);
5173
237
    if (add_lexmin_eq(tab, v->el) < 0)
5174
0
      goto error;
5175
237
  }
5176
123
5177
123
  isl_vec_free(v);
5178
123
  return n;
5179
0
error:
5180
0
  isl_vec_free(v);
5181
0
  return -1;
5182
123
}
5183
5184
/* Global internal data for isl_tab_basic_set_non_trivial_lexmin.
5185
 *
5186
 * "v" is a pre-allocated vector that can be used for adding
5187
 * constraints to the tableau.
5188
 */
5189
struct isl_trivial_global {
5190
  isl_vec *v;
5191
};
5192
5193
/* Fix triviality direction "dir" of the given region to zero.
5194
 *
5195
 * This function assumes that at least two more rows and at least
5196
 * two more elements in the constraint array are available in the tableau.
5197
 */
5198
static isl_stat fix_zero(struct isl_tab *tab, struct isl_trivial_region *region,
5199
  int dir, struct isl_trivial_global *data)
5200
64
{
5201
64
  int len;
5202
64
5203
64
  data->v = isl_vec_clr(data->v);
5204
64
  if (!data->v)
5205
0
    return isl_stat_error;
5206
64
  len = isl_mat_cols(region->trivial);
5207
64
  isl_seq_cpy(data->v->el + 1 + region->pos, region->trivial->row[dir],
5208
64
        len);
5209
64
  if (add_lexmin_eq(tab, data->v->el) < 0)
5210
0
    return isl_stat_error;
5211
64
5212
64
  return isl_stat_ok;
5213
64
}
5214
5215
/* This function selects case "side" for non-triviality region "region",
5216
 * assuming all the equality constraints have been imposed already.
5217
 * In particular, the triviality direction side/2 is made positive
5218
 * if side is even and made negative if side is odd.
5219
 *
5220
 * This function assumes that at least one more row and at least
5221
 * one more element in the constraint array are available in the tableau.
5222
 */
5223
static struct isl_tab *pos_neg(struct isl_tab *tab,
5224
  struct isl_trivial_region *region,
5225
  int side, struct isl_trivial_global *data)
5226
783
{
5227
783
  int len;
5228
783
5229
783
  data->v = isl_vec_clr(data->v);
5230
783
  if (!data->v)
5231
0
    goto error;
5232
783
  
isl_int_set_si783
(data->v->el[0], -1);
5233
783
  len = isl_mat_cols(region->trivial);
5234
783
  if (side % 2 == 0)
5235
487
    isl_seq_cpy(data->v->el + 1 + region->pos,
5236
487
          region->trivial->row[side / 2], len);
5237
783
  else
5238
296
    isl_seq_neg(data->v->el + 1 + region->pos,
5239
296
          region->trivial->row[side / 2], len);
5240
783
  return add_lexmin_ineq(tab, data->v->el);
5241
0
error:
5242
0
  isl_tab_free(tab);
5243
0
  return NULL;
5244
783
}
5245
5246
/* Local data at each level of the backtracking procedure of
5247
 * isl_tab_basic_set_non_trivial_lexmin.
5248
 *
5249
 * "update" is set if a solution has been found in the current case
5250
 * of this level, such that a better solution needs to be enforced
5251
 * in the next case.
5252
 * "n_zero" is the number of initial coordinates that have already
5253
 * been forced to be zero at this level.
5254
 * "region" is the non-triviality region considered at this level.
5255
 * "side" is the index of the current case at this level.
5256
 * "n" is the number of triviality directions.
5257
 */
5258
struct isl_trivial {
5259
  int update;
5260
  int n_zero;
5261
  int region;
5262
  int side;
5263
  int n;
5264
  struct isl_tab_undo *snap;
5265
};
5266
5267
/* Return the lexicographically smallest non-trivial solution of the
5268
 * given ILP problem.
5269
 *
5270
 * All variables are assumed to be non-negative.
5271
 *
5272
 * n_op is the number of initial coordinates to optimize.
5273
 * That is, once a solution has been found, we will only continue looking
5274
 * for solutions that result in significantly better values for those
5275
 * initial coordinates.  That is, we only continue looking for solutions
5276
 * that increase the number of initial zeros in this sequence.
5277
 *
5278
 * A solution is non-trivial, if it is non-trivial on each of the
5279
 * specified regions.  Each region represents a sequence of
5280
 * triviality directions on a sequence of variables that starts
5281
 * at a given position.  A solution is non-trivial on such a region if
5282
 * at least one of the triviality directions is non-zero
5283
 * on that sequence of variables.
5284
 *
5285
 * Whenever a conflict is encountered, all constraints involved are
5286
 * reported to the caller through a call to "conflict".
5287
 *
5288
 * We perform a simple branch-and-bound backtracking search.
5289
 * Each level in the search represents an initially trivial region
5290
 * that is forced to be non-trivial.
5291
 * At each level we consider 2 * n cases, where n
5292
 * is the number of triviality directions.
5293
 * In terms of those n directions v_i, we consider the cases
5294
 *  v_0 >= 1
5295
 *  v_0 <= -1
5296
 *  v_0 = 0 and v_1 >= 1
5297
 *  v_0 = 0 and v_1 <= -1
5298
 *  v_0 = 0 and v_1 = 0 and v_2 >= 1
5299
 *  v_0 = 0 and v_1 = 0 and v_2 <= -1
5300
 *  ...
5301
 * in this order.
5302
 */
5303
__isl_give isl_vec *isl_tab_basic_set_non_trivial_lexmin(
5304
  __isl_take isl_basic_set *bset, int n_op, int n_region,
5305
  struct isl_trivial_region *region,
5306
  int (*conflict)(int con, void *user), void *user)
5307
403
{
5308
403
  struct isl_trivial_global data = { 0 };
5309
403
  int i;
5310
403
  int r;
5311
403
  isl_ctx *ctx;
5312
403
  isl_vec *sol = NULL;
5313
403
  struct isl_tab *tab;
5314
403
  struct isl_trivial *triv = NULL;
5315
403
  int level, init;
5316
403
5317
403
  if (!bset)
5318
0
    return NULL;
5319
403
5320
403
  ctx = isl_basic_set_get_ctx(bset);
5321
403
  sol = isl_vec_alloc(ctx, 0);
5322
403
5323
403
  tab = tab_for_lexmin(bset, NULL, 0, 0);
5324
403
  if (!tab)
5325
0
    goto error;
5326
403
  tab->conflict = conflict;
5327
403
  tab->conflict_user = user;
5328
403
5329
403
  data.v = isl_vec_alloc(ctx, 1 + tab->n_var);
5330
403
  triv = isl_calloc_array(ctx, struct isl_trivial, n_region);
5331
403
  if (
!data.v || 403
(n_region && 403
!triv403
))
5332
0
    goto error;
5333
403
5334
403
  level = 0;
5335
403
  init = 1;
5336
403
5337
1.99k
  while (
level >= 01.99k
) {
5338
1.77k
    int side, base;
5339
1.77k
5340
1.77k
    if (
init1.77k
) {
5341
1.18k
      tab = cut_to_integer_lexmin(tab, CUT_ONE);
5342
1.18k
      if (!tab)
5343
0
        goto error;
5344
1.18k
      
if (1.18k
tab->empty1.18k
)
5345
465
        goto backtrack;
5346
721
      r = first_trivial_region(tab, n_region, region);
5347
721
      if (r < 0)
5348
0
        goto error;
5349
721
      
if (721
r == n_region721
) {
5350
613
        for (i = 0; 
i < level613
;
++i315
)
5351
315
          triv[i].update = 1;
5352
298
        isl_vec_free(sol);
5353
298
        sol = isl_tab_get_sample_value(tab);
5354
298
        if (!sol)
5355
0
          goto error;
5356
298
        
if (298
is_optimal(sol, n_op)298
)
5357
186
          break;
5358
112
        goto backtrack;
5359
112
      }
5360
423
      
if (423
level >= n_region423
)
5361
0
        isl_die(ctx, isl_error_internal,
5362
423
          "nesting level too deep", goto error);
5363
423
      triv[level].n = isl_mat_rows(region[r].trivial);
5364
423
      if (isl_tab_extend_cons(tab,
5365
423
              2 * triv[level].n + 2 * n_op) < 0)
5366
0
        goto error;
5367
423
      triv[level].region = r;
5368
423
      triv[level].side = 0;
5369
423
      triv[level].update = 0;
5370
423
      triv[level].n_zero = 0;
5371
423
    }
5372
1.77k
5373
1.01k
    r = triv[level].region;
5374
1.01k
    side = triv[level].side;
5375
1.01k
    base = 2 * (side/2);
5376
1.01k
5377
1.01k
    if (
side >= 2 * triv[level].n1.01k
) {
5378
808
backtrack:
5379
808
      level--;
5380
808
      init = 0;
5381
808
      if (level >= 0)
5382
591
        
if (591
isl_tab_rollback(tab, triv[level].snap) < 0591
)
5383
0
          goto error;
5384
808
      continue;
5385
808
    }
5386
783
5387
783
    
if (783
triv[level].update783
) {
5388
123
      triv[level].n_zero = force_better_solution(tab, sol,
5389
123
                n_op, triv[level].n_zero);
5390
123
      if (triv[level].n_zero < 0)
5391
0
        goto error;
5392
123
      triv[level].update = 0;
5393
123
    }
5394
783
5395
783
    
if (783
side == base && 783
base >= 2487
&&
5396
64
        fix_zero(tab, &region[r], base / 2 - 1, &data) < 0)
5397
0
      goto error;
5398
783
5399
783
    triv[level].snap = isl_tab_snap(tab);
5400
783
    if (isl_tab_push_basis(tab) < 0)
5401
0
      goto error;
5402
783
5403
783
    tab = pos_neg(tab, &region[r], side, &data);
5404
783
    if (!tab)
5405
0
      goto error;
5406
783
5407
783
    triv[level].side++;
5408
783
    level++;
5409
783
    init = 1;
5410
783
  }
5411
403
5412
403
  free(triv);
5413
403
  isl_vec_free(data.v);
5414
403
  isl_tab_free(tab);
5415
403
  isl_basic_set_free(bset);
5416
403
5417
403
  return sol;
5418
0
error:
5419
0
  free(triv);
5420
0
  isl_vec_free(data.v);
5421
0
  isl_tab_free(tab);
5422
0
  isl_basic_set_free(bset);
5423
0
  isl_vec_free(sol);
5424
0
  return NULL;
5425
403
}
5426
5427
/* Wrapper for a tableau that is used for computing
5428
 * the lexicographically smallest rational point of a non-negative set.
5429
 * This point is represented by the sample value of "tab",
5430
 * unless "tab" is empty.
5431
 */
5432
struct isl_tab_lexmin {
5433
  isl_ctx *ctx;
5434
  struct isl_tab *tab;
5435
};
5436
5437
/* Free "tl" and return NULL.
5438
 */
5439
__isl_null isl_tab_lexmin *isl_tab_lexmin_free(__isl_take isl_tab_lexmin *tl)
5440
32
{
5441
32
  if (!tl)
5442
0
    return NULL;
5443
32
  isl_ctx_deref(tl->ctx);
5444
32
  isl_tab_free(tl->tab);
5445
32
  free(tl);
5446
32
5447
32
  return NULL;
5448
32
}
5449
5450
/* Construct an isl_tab_lexmin for computing
5451
 * the lexicographically smallest rational point in "bset",
5452
 * assuming that all variables are non-negative.
5453
 */
5454
__isl_give isl_tab_lexmin *isl_tab_lexmin_from_basic_set(
5455
  __isl_take isl_basic_set *bset)
5456
32
{
5457
32
  isl_ctx *ctx;
5458
32
  isl_tab_lexmin *tl;
5459
32
5460
32
  if (!bset)
5461
0
    return NULL;
5462
32
5463
32
  ctx = isl_basic_set_get_ctx(bset);
5464
32
  tl = isl_calloc_type(ctx, struct isl_tab_lexmin);
5465
32
  if (!tl)
5466
0
    goto error;
5467
32
  tl->ctx = ctx;
5468
32
  isl_ctx_ref(ctx);
5469
32
  tl->tab = tab_for_lexmin(bset, NULL, 0, 0);
5470
32
  isl_basic_set_free(bset);
5471
32
  if (!tl->tab)
5472
0
    return isl_tab_lexmin_free(tl);
5473
32
  return tl;
5474
0
error:
5475
0
  isl_basic_set_free(bset);
5476
0
  isl_tab_lexmin_free(tl);
5477
0
  return NULL;
5478
32
}
5479
5480
/* Return the dimension of the set represented by "tl".
5481
 */
5482
int isl_tab_lexmin_dim(__isl_keep isl_tab_lexmin *tl)
5483
4
{
5484
4
  return tl ? 
tl->tab->n_var4
:
-10
;
5485
4
}
5486
5487
/* Add the equality with coefficients "eq" to "tl", updating the optimal
5488
 * solution if needed.
5489
 * The equality is added as two opposite inequality constraints.
5490
 */
5491
__isl_give isl_tab_lexmin *isl_tab_lexmin_add_eq(__isl_take isl_tab_lexmin *tl,
5492
  isl_int *eq)
5493
4
{
5494
4
  unsigned n_var;
5495
4
5496
4
  if (
!tl || 4
!eq4
)
5497
0
    return isl_tab_lexmin_free(tl);
5498
4
5499
4
  
if (4
isl_tab_extend_cons(tl->tab, 2) < 04
)
5500
0
    return isl_tab_lexmin_free(tl);
5501
4
  n_var = tl->tab->n_var;
5502
4
  isl_seq_neg(eq, eq, 1 + n_var);
5503
4
  tl->tab = add_lexmin_ineq(tl->tab, eq);
5504
4
  isl_seq_neg(eq, eq, 1 + n_var);
5505
4
  tl->tab = add_lexmin_ineq(tl->tab, eq);
5506
4
5507
4
  if (!tl->tab)
5508
0
    return isl_tab_lexmin_free(tl);
5509
4
5510
4
  return tl;
5511
4
}
5512
5513
/* Add cuts to "tl" until the sample value reaches an integer value or
5514
 * until the result becomes empty.
5515
 */
5516
__isl_give isl_tab_lexmin *isl_tab_lexmin_cut_to_integer(
5517
  __isl_take isl_tab_lexmin *tl)
5518
2
{
5519
2
  if (!tl)
5520
0
    return NULL;
5521
2
  
tl->tab = cut_to_integer_lexmin(tl->tab, 2
CUT_ONE2
);
5522
2
  if (!tl->tab)
5523
0
    return isl_tab_lexmin_free(tl);
5524
2
  return tl;
5525
2
}
5526
5527
/* Return the lexicographically smallest rational point in the basic set
5528
 * from which "tl" was constructed.
5529
 * If the original input was empty, then return a zero-length vector.
5530
 */
5531
__isl_give isl_vec *isl_tab_lexmin_get_solution(__isl_keep isl_tab_lexmin *tl)
5532
38
{
5533
38
  if (!tl)
5534
0
    return NULL;
5535
38
  
if (38
tl->tab->empty38
)
5536
0
    return isl_vec_alloc(tl->ctx, 0);
5537
38
  else
5538
38
    return isl_tab_get_sample_value(tl->tab);
5539
0
}
5540
5541
struct isl_sol_pma {
5542
  struct isl_sol  sol;
5543
  isl_pw_multi_aff *pma;
5544
  isl_set *empty;
5545
};
5546
5547
static void sol_pma_free(struct isl_sol *sol)
5548
3.45k
{
5549
3.45k
  struct isl_sol_pma *sol_pma = (struct isl_sol_pma *) sol;
5550
3.45k
  isl_pw_multi_aff_free(sol_pma->pma);
5551
3.45k
  isl_set_free(sol_pma->empty);
5552
3.45k
}
5553
5554
/* This function is called for parts of the context where there is
5555
 * no solution, with "bset" corresponding to the context tableau.
5556
 * Simply add the basic set to the set "empty".
5557
 */
5558
static void sol_pma_add_empty(struct isl_sol_pma *sol,
5559
  __isl_take isl_basic_set *bset)
5560
1.25k
{
5561
1.25k
  if (
!bset || 1.25k
!sol->empty1.25k
)
5562
0
    goto error;
5563
1.25k
5564
1.25k
  sol->empty = isl_set_grow(sol->empty, 1);
5565
1.25k
  bset = isl_basic_set_simplify(bset);
5566
1.25k
  bset = isl_basic_set_finalize(bset);
5567
1.25k
  sol->empty = isl_set_add_basic_set(sol->empty, bset);
5568
1.25k
  if (!sol->empty)
5569
0
    sol->sol.error = 1;
5570
1.25k
  return;
5571
0
error:
5572
0
  isl_basic_set_free(bset);
5573
0
  sol->sol.error = 1;
5574
0
}
5575
5576
/* Given a basic set "dom" that represents the context and a tuple of
5577
 * affine expressions "maff" defined over this domain, construct
5578
 * an isl_pw_multi_aff with a single cell corresponding to "dom" and
5579
 * the affine expressions in "maff".
5580
 */
5581
static void sol_pma_add(struct isl_sol_pma *sol,
5582
  __isl_take isl_basic_set *dom, __isl_take isl_multi_aff *maff)
5583
2.71k
{
5584
2.71k
  isl_pw_multi_aff *pma;
5585
2.71k
5586
2.71k
  dom = isl_basic_set_simplify(dom);
5587
2.71k
  dom = isl_basic_set_finalize(dom);
5588
2.71k
  pma = isl_pw_multi_aff_alloc(isl_set_from_basic_set(dom), maff);
5589
2.71k
  sol->pma = isl_pw_multi_aff_add_disjoint(sol->pma, pma);
5590
2.71k
  if (!sol->pma)
5591
0
    sol->sol.error = 1;
5592
2.71k
}
5593
5594
static void sol_pma_add_empty_wrap(struct isl_sol *sol,
5595
  __isl_take isl_basic_set *bset)
5596
1.25k
{
5597
1.25k
  sol_pma_add_empty((struct isl_sol_pma *)sol, bset);
5598
1.25k
}
5599
5600
static void sol_pma_add_wrap(struct isl_sol *sol,
5601
  __isl_take isl_basic_set *dom, __isl_take isl_multi_aff *ma)
5602
2.71k
{
5603
2.71k
  sol_pma_add((struct isl_sol_pma *)sol, dom, ma);
5604
2.71k
}
5605
5606
/* Construct an isl_sol_pma structure for accumulating the solution.
5607
 * If track_empty is set, then we also keep track of the parts
5608
 * of the context where there is no solution.
5609
 * If max is set, then we are solving a maximization, rather than
5610
 * a minimization problem, which means that the variables in the
5611
 * tableau have value "M - x" rather than "M + x".
5612
 */
5613
static struct isl_sol *sol_pma_init(__isl_keep isl_basic_map *bmap,
5614
  __isl_take isl_basic_set *dom, int track_empty, int max)
5615
3.45k
{
5616
3.45k
  struct isl_sol_pma *sol_pma = NULL;
5617
3.45k
  isl_space *space;
5618
3.45k
5619
3.45k
  if (!bmap)
5620
0
    goto error;
5621
3.45k
5622
3.45k
  
sol_pma = 3.45k
isl_calloc_type3.45k
(bmap->ctx, struct isl_sol_pma);
5623
3.45k
  if (!sol_pma)
5624
0
    goto error;
5625
3.45k
5626
3.45k
  sol_pma->sol.free = &sol_pma_free;
5627
3.45k
  if (sol_init(&sol_pma->sol, bmap, dom, max) < 0)
5628
0
    goto error;
5629
3.45k
  sol_pma->sol.add = &sol_pma_add_wrap;
5630
1.08k
  sol_pma->sol.add_empty = track_empty ? &sol_pma_add_empty_wrap : NULL;
5631
3.45k
  space = isl_space_copy(sol_pma->sol.space);
5632
3.45k
  sol_pma->pma = isl_pw_multi_aff_empty(space);
5633
3.45k
  if (!sol_pma->pma)
5634
0
    goto error;
5635
3.45k
5636
3.45k
  
if (3.45k
track_empty3.45k
) {
5637
1.08k
    sol_pma->empty = isl_set_alloc_space(isl_basic_set_get_space(dom),
5638
1.08k
              1, ISL_SET_DISJOINT);
5639
1.08k
    if (!sol_pma->empty)
5640
0
      goto error;
5641
3.45k
  }
5642
3.45k
5643
3.45k
  isl_basic_set_free(dom);
5644
3.45k
  return &sol_pma->sol;
5645
0
error:
5646
0
  isl_basic_set_free(dom);
5647
0
  sol_free(&sol_pma->sol);
5648
0
  return NULL;
5649
3.45k
}
5650
5651
/* Base case of isl_tab_basic_map_partial_lexopt, after removing
5652
 * some obvious symmetries.
5653
 *
5654
 * We call basic_map_partial_lexopt_base_sol and extract the results.
5655
 */
5656
static __isl_give isl_pw_multi_aff *basic_map_partial_lexopt_base_pw_multi_aff(
5657
  __isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom,
5658
  __isl_give isl_set **empty, int max)
5659
3.45k
{
5660
3.45k
  isl_pw_multi_aff *result = NULL;
5661
3.45k
  struct isl_sol *sol;
5662
3.45k
  struct isl_sol_pma *sol_pma;
5663
3.45k
5664
3.45k
  sol = basic_map_partial_lexopt_base_sol(bmap, dom, empty, max,
5665
3.45k
            &sol_pma_init);
5666
3.45k
  if (!sol)
5667
0
    return NULL;
5668
3.45k
  sol_pma = (struct isl_sol_pma *) sol;
5669
3.45k
5670
3.45k
  result = isl_pw_multi_aff_copy(sol_pma->pma);
5671
3.45k
  if (empty)
5672
1.08k
    *empty = isl_set_copy(sol_pma->empty);
5673
3.45k
  sol_free(&sol_pma->sol);
5674
3.45k
  return result;
5675
3.45k
}
5676
5677
/* Given that the last input variable of "maff" represents the minimum
5678
 * of some bounds, check whether we need to plug in the expression
5679
 * of the minimum.
5680
 *
5681
 * In particular, check if the last input variable appears in any
5682
 * of the expressions in "maff".
5683
 */
5684
static int need_substitution(__isl_keep isl_multi_aff *maff)
5685
113
{
5686
113
  int i;
5687
113
  unsigned pos;
5688
113
5689
113
  pos = isl_multi_aff_dim(maff, isl_dim_in) - 1;
5690
113
5691
166
  for (i = 0; 
i < maff->n166
;
++i53
)
5692
159
    
if (159
isl_aff_involves_dims(maff->p[i], isl_dim_in, pos, 1)159
)
5693
106
      return 1;
5694
113
5695
7
  return 0;
5696
113
}
5697
5698
/* Given a set of upper bounds on the last "input" variable m,
5699
 * construct a piecewise affine expression that selects
5700
 * the minimal upper bound to m, i.e.,
5701
 * divide the space into cells where one
5702
 * of the upper bounds is smaller than all the others and select
5703
 * this upper bound on that cell.
5704
 *
5705
 * In particular, if there are n bounds b_i, then the result
5706
 * consists of n cell, each one of the form
5707
 *
5708
 *  b_i <= b_j  for j > i
5709
 *  b_i <  b_j  for j < i
5710
 *
5711
 * The affine expression on this cell is
5712
 *
5713
 *  b_i
5714
 */
5715
static __isl_give isl_pw_aff *set_minimum_pa(__isl_take isl_space *space,
5716
  __isl_take isl_mat *var)
5717
110
{
5718
110
  int i;
5719
110
  isl_aff *aff = NULL;
5720
110
  isl_basic_set *bset = NULL;
5721
110
  isl_pw_aff *paff = NULL;
5722
110
  isl_space *pw_space;
5723
110
  isl_local_space *ls = NULL;
5724
110
5725
110
  if (
!space || 110
!var110
)
5726
0
    goto error;
5727
110
5728
110
  ls = isl_local_space_from_space(isl_space_copy(space));
5729
110
  pw_space = isl_space_copy(space);
5730
110
  pw_space = isl_space_from_domain(pw_space);
5731
110
  pw_space = isl_space_add_dims(pw_space, isl_dim_out, 1);
5732
110
  paff = isl_pw_aff_alloc_size(pw_space, var->n_row);
5733
110
5734
333
  for (i = 0; 
i < var->n_row333
;
++i223
) {
5735
223
    isl_pw_aff *paff_i;
5736
223
5737
223
    aff = isl_aff_alloc(isl_local_space_copy(ls));
5738
223
    bset = isl_basic_set_alloc_space(isl_space_copy(space), 0,
5739
223
                 0, var->n_row - 1);
5740
223
    if (
!aff || 223
!bset223
)
5741
0
      goto error;
5742
223
    
isl_int_set_si223
(aff->v->el[0], 1);
5743
223
    isl_seq_cpy(aff->v->el + 1, var->row[i], var->n_col);
5744
223
    isl_int_set_si(aff->v->el[1 + var->n_col], 0);
5745
223
    bset = select_minimum(bset, var, i);
5746
223
    paff_i = isl_pw_aff_alloc(isl_set_from_basic_set(bset), aff);
5747
223
    paff = isl_pw_aff_add_disjoint(paff, paff_i);
5748
223
  }
5749
110
5750
110
  isl_local_space_free(ls);
5751
110
  isl_space_free(space);
5752
110
  isl_mat_free(var);
5753
110
  return paff;
5754
0
error:
5755
0
  isl_aff_free(aff);
5756
0
  isl_basic_set_free(bset);
5757
0
  isl_pw_aff_free(paff);
5758
0
  isl_local_space_free(ls);
5759
0
  isl_space_free(space);
5760
0
  isl_mat_free(var);
5761
0
  return NULL;
5762
110
}
5763
5764
/* Given a piecewise multi-affine expression of which the last input variable
5765
 * is the minimum of the bounds in "cst", plug in the value of the minimum.
5766
 * This minimum expression is given in "min_expr_pa".
5767
 * The set "min_expr" contains the same information, but in the form of a set.
5768
 * The variable is subsequently projected out.
5769
 *
5770
 * The implementation is similar to those of "split" and "split_domain".
5771
 * If the variable appears in a given expression, then minimum expression
5772
 * is plugged in.  Otherwise, if the variable appears in the constraints
5773
 * and a split is required, then the domain is split.  Otherwise, no split
5774
 * is performed.
5775
 */
5776
static __isl_give isl_pw_multi_aff *split_domain_pma(
5777
  __isl_take isl_pw_multi_aff *opt, __isl_take isl_pw_aff *min_expr_pa,
5778
  __isl_take isl_set *min_expr, __isl_take isl_mat *cst)
5779
110
{
5780
110
  int n_in;
5781
110
  int i;
5782
110
  isl_space *space;
5783
110
  isl_pw_multi_aff *res;
5784
110
5785
110
  if (
!opt || 110
!min_expr110
||
!cst110
)
5786
0
    goto error;
5787
110
5788
110
  n_in = isl_pw_multi_aff_dim(opt, isl_dim_in);
5789
110
  space = isl_pw_multi_aff_get_space(opt);
5790
110
  space = isl_space_drop_dims(space, isl_dim_in, n_in - 1, 1);
5791
110
  res = isl_pw_multi_aff_empty(space);
5792
110
5793
223
  for (i = 0; 
i < opt->n223
;
++i113
) {
5794
113
    isl_pw_multi_aff *pma;
5795
113
5796
113
    pma = isl_pw_multi_aff_alloc(isl_set_copy(opt->p[i].set),
5797
113
           isl_multi_aff_copy(opt->p[i].maff));
5798
113
    if (need_substitution(opt->p[i].maff))
5799
106
      pma = isl_pw_multi_aff_substitute(pma,
5800
106
          isl_dim_in, n_in - 1, min_expr_pa);
5801
7
    else {
5802
7
      isl_bool split;
5803
7
      split = need_split_set(opt->p[i].set, cst);
5804
7
      if (split < 0)
5805
0
        pma = isl_pw_multi_aff_free(pma);
5806
7
      else 
if (7
split7
)
5807
0
        pma = isl_pw_multi_aff_intersect_domain(pma,
5808
0
                   isl_set_copy(min_expr));
5809
7
    }
5810
113
    pma = isl_pw_multi_aff_project_out(pma,
5811
113
                isl_dim_in, n_in - 1, 1);
5812
113
5813
113
    res = isl_pw_multi_aff_add_disjoint(res, pma);
5814
113
  }
5815
110
5816
110
  isl_pw_multi_aff_free(opt);
5817
110
  isl_pw_aff_free(min_expr_pa);
5818
110
  isl_set_free(min_expr);
5819
110
  isl_mat_free(cst);
5820
110
  return res;
5821
0
error:
5822
0
  isl_pw_multi_aff_free(opt);
5823
0
  isl_pw_aff_free(min_expr_pa);
5824
0
  isl_set_free(min_expr);
5825
0
  isl_mat_free(cst);
5826
0
  return NULL;
5827
110
}
5828
5829
static __isl_give isl_pw_multi_aff *basic_map_partial_lexopt_pw_multi_aff(
5830
  __isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom,
5831
  __isl_give isl_set **empty, int max);
5832
5833
/* This function is called from basic_map_partial_lexopt_symm.
5834
 * The last variable of "bmap" and "dom" corresponds to the minimum
5835
 * of the bounds in "cst".  "map_space" is the space of the original
5836
 * input relation (of basic_map_partial_lexopt_symm) and "set_space"
5837
 * is the space of the original domain.
5838
 *
5839
 * We recursively call basic_map_partial_lexopt and then plug in
5840
 * the definition of the minimum in the result.
5841
 */
5842
static __isl_give isl_pw_multi_aff *
5843
basic_map_partial_lexopt_symm_core_pw_multi_aff(
5844
  __isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom,
5845
  __isl_give isl_set **empty, int max, __isl_take isl_mat *cst,
5846
  __isl_take isl_space *map_space, __isl_take isl_space *set_space)
5847
110
{
5848
110
  isl_pw_multi_aff *opt;
5849
110
  isl_pw_aff *min_expr_pa;
5850
110
  isl_set *min_expr;
5851
110
5852
110
  min_expr = set_minimum(isl_basic_set_get_space(dom), isl_mat_copy(cst));
5853
110
  min_expr_pa = set_minimum_pa(isl_basic_set_get_space(dom),
5854
110
          isl_mat_copy(cst));
5855
110
5856
110
  opt = basic_map_partial_lexopt_pw_multi_aff(bmap, dom, empty, max);
5857
110
5858
110
  if (
empty110
) {
5859
0
    *empty = split(*empty,
5860
0
             isl_set_copy(min_expr), isl_mat_copy(cst));
5861
0
    *empty = isl_set_reset_space(*empty, set_space);
5862
0
  }
5863
110
5864
110
  opt = split_domain_pma(opt, min_expr_pa, min_expr, cst);
5865
110
  opt = isl_pw_multi_aff_reset_space(opt, map_space);
5866
110
5867
110
  return opt;
5868
110
}
5869
5870
#undef TYPE
5871
#define TYPE  isl_pw_multi_aff
5872
#undef SUFFIX
5873
#define SUFFIX  _pw_multi_aff
5874
#include "isl_tab_lexopt_templ.c"