Coverage Report

Created: 2017-11-23 03:11

/Users/buildslave/jenkins/workspace/clang-stage2-coverage-R/llvm/tools/polly/lib/External/isl/isl_tab_pip.c
Line
Count
Source (jump to first uncovered line)
1
/*
2
 * Copyright 2008-2009 Katholieke Universiteit Leuven
3
 * Copyright 2010      INRIA Saclay
4
 * Copyright 2016-2017 Sven Verdoolaege
5
 *
6
 * Use of this software is governed by the MIT license
7
 *
8
 * Written by Sven Verdoolaege, K.U.Leuven, Departement
9
 * Computerwetenschappen, Celestijnenlaan 200A, B-3001 Leuven, Belgium
10
 * and INRIA Saclay - Ile-de-France, Parc Club Orsay Universite,
11
 * ZAC des vignes, 4 rue Jacques Monod, 91893 Orsay, France 
12
 */
13
14
#include <isl_ctx_private.h>
15
#include "isl_map_private.h"
16
#include <isl_seq.h>
17
#include "isl_tab.h"
18
#include "isl_sample.h"
19
#include <isl_mat_private.h>
20
#include <isl_vec_private.h>
21
#include <isl_aff_private.h>
22
#include <isl_constraint_private.h>
23
#include <isl_options_private.h>
24
#include <isl_config.h>
25
26
#include <bset_to_bmap.c>
27
28
/*
29
 * The implementation of parametric integer linear programming in this file
30
 * was inspired by the paper "Parametric Integer Programming" and the
31
 * report "Solving systems of affine (in)equalities" by Paul Feautrier
32
 * (and others).
33
 *
34
 * The strategy used for obtaining a feasible solution is different
35
 * from the one used in isl_tab.c.  In particular, in isl_tab.c,
36
 * upon finding a constraint that is not yet satisfied, we pivot
37
 * in a row that increases the constant term of the row holding the
38
 * constraint, making sure the sample solution remains feasible
39
 * for all the constraints it already satisfied.
40
 * Here, we always pivot in the row holding the constraint,
41
 * choosing a column that induces the lexicographically smallest
42
 * increment to the sample solution.
43
 *
44
 * By starting out from a sample value that is lexicographically
45
 * smaller than any integer point in the problem space, the first
46
 * feasible integer sample point we find will also be the lexicographically
47
 * smallest.  If all variables can be assumed to be non-negative,
48
 * then the initial sample value may be chosen equal to zero.
49
 * However, we will not make this assumption.  Instead, we apply
50
 * the "big parameter" trick.  Any variable x is then not directly
51
 * used in the tableau, but instead it is represented by another
52
 * variable x' = M + x, where M is an arbitrarily large (positive)
53
 * value.  x' is therefore always non-negative, whatever the value of x.
54
 * Taking as initial sample value x' = 0 corresponds to x = -M,
55
 * which is always smaller than any possible value of x.
56
 *
57
 * The big parameter trick is used in the main tableau and
58
 * also in the context tableau if isl_context_lex is used.
59
 * In this case, each tableaus has its own big parameter.
60
 * Before doing any real work, we check if all the parameters
61
 * happen to be non-negative.  If so, we drop the column corresponding
62
 * to M from the initial context tableau.
63
 * If isl_context_gbr is used, then the big parameter trick is only
64
 * used in the main tableau.
65
 */
66
67
struct isl_context;
68
struct isl_context_op {
69
  /* detect nonnegative parameters in context and mark them in tab */
70
  struct isl_tab *(*detect_nonnegative_parameters)(
71
      struct isl_context *context, struct isl_tab *tab);
72
  /* return temporary reference to basic set representation of context */
73
  struct isl_basic_set *(*peek_basic_set)(struct isl_context *context);
74
  /* return temporary reference to tableau representation of context */
75
  struct isl_tab *(*peek_tab)(struct isl_context *context);
76
  /* add equality; check is 1 if eq may not be valid;
77
   * update is 1 if we may want to call ineq_sign on context later.
78
   */
79
  void (*add_eq)(struct isl_context *context, isl_int *eq,
80
      int check, int update);
81
  /* add inequality; check is 1 if ineq may not be valid;
82
   * update is 1 if we may want to call ineq_sign on context later.
83
   */
84
  void (*add_ineq)(struct isl_context *context, isl_int *ineq,
85
      int check, int update);
86
  /* check sign of ineq based on previous information.
87
   * strict is 1 if saturation should be treated as a positive sign.
88
   */
89
  enum isl_tab_row_sign (*ineq_sign)(struct isl_context *context,
90
      isl_int *ineq, int strict);
91
  /* check if inequality maintains feasibility */
92
  int (*test_ineq)(struct isl_context *context, isl_int *ineq);
93
  /* return index of a div that corresponds to "div" */
94
  int (*get_div)(struct isl_context *context, struct isl_tab *tab,
95
      struct isl_vec *div);
96
  /* insert div "div" to context at "pos" and return non-negativity */
97
  isl_bool (*insert_div)(struct isl_context *context, int pos,
98
    __isl_keep isl_vec *div);
99
  int (*detect_equalities)(struct isl_context *context,
100
      struct isl_tab *tab);
101
  /* return row index of "best" split */
102
  int (*best_split)(struct isl_context *context, struct isl_tab *tab);
103
  /* check if context has already been determined to be empty */
104
  int (*is_empty)(struct isl_context *context);
105
  /* check if context is still usable */
106
  int (*is_ok)(struct isl_context *context);
107
  /* save a copy/snapshot of context */
108
  void *(*save)(struct isl_context *context);
109
  /* restore saved context */
110
  void (*restore)(struct isl_context *context, void *);
111
  /* discard saved context */
112
  void (*discard)(void *);
113
  /* invalidate context */
114
  void (*invalidate)(struct isl_context *context);
115
  /* free context */
116
  __isl_null struct isl_context *(*free)(struct isl_context *context);
117
};
118
119
/* Shared parts of context representation.
120
 *
121
 * "n_unknown" is the number of final unknown integer divisions
122
 * in the input domain.
123
 */
124
struct isl_context {
125
  struct isl_context_op *op;
126
  int n_unknown;
127
};
128
129
struct isl_context_lex {
130
  struct isl_context context;
131
  struct isl_tab *tab;
132
};
133
134
/* A stack (linked list) of solutions of subtrees of the search space.
135
 *
136
 * "ma" describes the solution as a function of "dom".
137
 * In particular, the domain space of "ma" is equal to the space of "dom".
138
 *
139
 * If "ma" is NULL, then there is no solution on "dom".
140
 */
141
struct isl_partial_sol {
142
  int level;
143
  struct isl_basic_set *dom;
144
  isl_multi_aff *ma;
145
146
  struct isl_partial_sol *next;
147
};
148
149
struct isl_sol;
150
struct isl_sol_callback {
151
  struct isl_tab_callback callback;
152
  struct isl_sol *sol;
153
};
154
155
/* isl_sol is an interface for constructing a solution to
156
 * a parametric integer linear programming problem.
157
 * Every time the algorithm reaches a state where a solution
158
 * can be read off from the tableau, the function "add" is called
159
 * on the isl_sol passed to find_solutions_main.  In a state where
160
 * the tableau is empty, "add_empty" is called instead.
161
 * "free" is called to free the implementation specific fields, if any.
162
 *
163
 * "error" is set if some error has occurred.  This flag invalidates
164
 * the remainder of the data structure.
165
 * If "rational" is set, then a rational optimization is being performed.
166
 * "level" is the current level in the tree with nodes for each
167
 * split in the context.
168
 * If "max" is set, then a maximization problem is being solved, rather than
169
 * a minimization problem, which means that the variables in the
170
 * tableau have value "M - x" rather than "M + x".
171
 * "n_out" is the number of output dimensions in the input.
172
 * "space" is the space in which the solution (and also the input) lives.
173
 *
174
 * The context tableau is owned by isl_sol and is updated incrementally.
175
 *
176
 * There are currently 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.98k
{
202
7.98k
  struct isl_partial_sol *partial, *next;
203
7.98k
  if (!sol)
204
0
    return;
205
7.98k
  for (partial = sol->partial; partial; 
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.98k
  isl_space_free(sol->space);
212
7.98k
  if (sol->context)
213
7.98k
    sol->context->op->free(sol->context);
214
7.98k
  sol->free(sol);
215
7.98k
  free(sol);
216
7.98k
}
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.58k
{
226
9.58k
  struct isl_partial_sol *partial;
227
9.58k
228
9.58k
  if (sol->error || !dom)
229
0
    goto error;
230
9.58k
231
9.58k
  partial = isl_alloc_type(dom->ctx, struct isl_partial_sol);
232
9.58k
  if (!partial)
233
0
    goto error;
234
9.58k
235
9.58k
  partial->level = sol->level;
236
9.58k
  partial->dom = dom;
237
9.58k
  partial->ma = ma;
238
9.58k
  partial->next = sol->partial;
239
9.58k
240
9.58k
  sol->partial = partial;
241
9.58k
242
9.58k
  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
7.16k
{
254
7.16k
  int i;
255
7.16k
  unsigned rows, cols, n;
256
7.16k
257
7.16k
  if (!M)
258
0
    return isl_stat_error;
259
7.16k
  rows = isl_mat_rows(M);
260
7.16k
  cols = isl_mat_cols(M);
261
7.16k
  n = cols - first;
262
30.6k
  for (i = 0; i < rows; 
++i23.4k
)
263
23.4k
    if (isl_seq_first_non_zero(M->row[i] + first, n) != -1)
264
23.4k
      
isl_die0
(isl_mat_get_ctx(M), isl_error_internal,
265
7.16k
        "final columns should be zero",
266
7.16k
        return isl_stat_error);
267
7.16k
  return isl_stat_ok;
268
7.16k
}
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
7.16k
{
279
7.16k
  int i, dim;
280
7.16k
  isl_aff *aff;
281
7.16k
282
7.16k
  if (!ma || !ls || !M)
283
0
    goto error;
284
7.16k
285
7.16k
  dim = isl_local_space_dim(ls, isl_dim_all);
286
7.16k
  if (check_final_columns_are_zero(M, 1 + dim) < 0)
287
0
    goto error;
288
23.4k
  
for (i = 1; 7.16k
i < M->n_row;
++i16.3k
) {
289
16.3k
    aff = isl_aff_alloc(isl_local_space_copy(ls));
290
16.3k
    if (aff) {
291
16.3k
      isl_int_set(aff->v->el[0], M->row[0][0]);
292
16.3k
      isl_seq_cpy(aff->v->el + 1, M->row[i], 1 + dim);
293
16.3k
    }
294
16.3k
    aff = isl_aff_normalize(aff);
295
16.3k
    ma = isl_multi_aff_set_aff(ma, i - 1, aff);
296
16.3k
  }
297
7.16k
  isl_local_space_free(ls);
298
7.16k
  isl_mat_free(M);
299
7.16k
300
7.16k
  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
7.16k
}
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
7.16k
{
331
7.16k
  isl_local_space *ls;
332
7.16k
  isl_multi_aff *ma;
333
7.16k
  int n_div, n_known;
334
7.16k
335
7.16k
  n_div = isl_basic_set_dim(dom, isl_dim_div);
336
7.16k
  n_known = n_div - sol->context->n_unknown;
337
7.16k
338
7.16k
  ma = isl_multi_aff_alloc(isl_space_copy(sol->space));
339
7.16k
  ls = isl_basic_set_get_local_space(dom);
340
7.16k
  ls = isl_local_space_drop_dims(ls, isl_dim_div,
341
7.16k
          n_known, n_div - n_known);
342
7.16k
  ma = set_from_affine_matrix(ma, ls, M);
343
7.16k
344
7.16k
  if (!ma)
345
0
    dom = isl_basic_set_free(dom);
346
7.16k
  sol_push_sol(sol, dom, ma);
347
7.16k
}
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.48k
{
354
9.48k
  struct isl_partial_sol *partial;
355
9.48k
356
9.48k
  partial = sol->partial;
357
9.48k
  sol->partial = partial->next;
358
9.48k
359
9.48k
  if (partial->ma)
360
7.06k
    sol->add(sol, partial->dom, partial->ma);
361
2.42k
  else
362
2.42k
    sol->add_empty(sol, partial->dom);
363
9.48k
  free(partial);
364
9.48k
}
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.68k
{
370
9.68k
  struct isl_basic_set *bset;
371
9.68k
372
9.68k
  if (sol->error)
373
0
    return NULL;
374
9.68k
375
9.68k
  bset = isl_basic_set_dup(sol->context->op->peek_basic_set(sol->context));
376
9.68k
  bset = isl_basic_set_update_from_tab(bset,
377
9.68k
      sol->context->op->peek_tab(sol->context));
378
9.68k
379
9.68k
  return bset;
380
9.68k
}
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.10k
{
387
2.10k
  if (!s1->ma != !s2->ma)
388
1.95k
    return isl_bool_false;
389
152
  if (!s1->ma)
390
0
    return isl_bool_true;
391
152
392
152
  return isl_multi_aff_plain_is_equal(s1->ma, s2->ma);
393
152
}
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
57
{
407
57
  struct isl_partial_sol *partial;
408
57
409
57
  partial = sol->partial;
410
57
  sol->partial = partial->next;
411
57
  partial->next = partial->next->next;
412
57
  sol->partial->next = partial;
413
57
}
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
99
{
427
99
  struct isl_partial_sol *partial;
428
99
  isl_basic_set *bset;
429
99
430
99
  partial = sol->partial;
431
99
432
99
  bset = sol_domain(sol);
433
99
  isl_basic_set_free(partial->next->dom);
434
99
  partial->next->dom = bset;
435
99
  partial->next->level = sol->level;
436
99
437
99
  if (!bset)
438
0
    return isl_stat_error;
439
99
440
99
  sol->partial = partial->next;
441
99
  isl_basic_set_free(partial->dom);
442
99
  isl_multi_aff_free(partial->ma);
443
99
  free(partial);
444
99
445
99
  return isl_stat_ok;
446
99
}
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
238
{
458
238
  isl_set *set;
459
238
  isl_pw_multi_aff *pma1, *pma2;
460
238
  isl_bool equal;
461
238
462
238
  set = isl_basic_set_compute_divs(isl_basic_set_copy(dom));
463
238
  pma1 = isl_pw_multi_aff_alloc(isl_set_copy(set),
464
238
          isl_multi_aff_copy(ma1));
465
238
  pma2 = isl_pw_multi_aff_alloc(set, isl_multi_aff_copy(ma2));
466
238
  equal = isl_pw_multi_aff_is_equal(pma1, pma2);
467
238
  isl_pw_multi_aff_free(pma1);
468
238
  isl_pw_multi_aff_free(pma2);
469
238
470
238
  return equal;
471
238
}
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.10k
{
490
2.10k
  struct isl_partial_sol *partial;
491
2.10k
  isl_bool same;
492
2.10k
493
2.10k
  partial = sol->partial;
494
2.10k
495
2.10k
  same = same_solution(partial, partial->next);
496
2.10k
  if (same < 0)
497
0
    return isl_stat_error;
498
2.10k
  if (same)
499
24
    return combine_initial_into_second(sol);
500
2.08k
  if (partial->ma && 
partial->next->ma146
) {
501
128
    same = equal_on_domain(partial->ma, partial->next->ma,
502
128
          partial->dom);
503
128
    if (same < 0)
504
0
      return isl_stat_error;
505
128
    if (same)
506
18
      return combine_initial_into_second(sol);
507
110
    same = equal_on_domain(partial->ma, partial->next->ma,
508
110
          partial->next->dom);
509
110
    if (same) {
510
57
      swap_initial(sol);
511
57
      return combine_initial_into_second(sol);
512
57
    }
513
2.00k
  }
514
2.00k
515
2.00k
  sol_pop_one(sol);
516
2.00k
  sol_pop_one(sol);
517
2.00k
518
2.00k
  return isl_stat_ok;
519
2.00k
}
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
9.36k
{
532
9.36k
  struct isl_partial_sol *partial;
533
9.36k
534
9.36k
  if (sol->error)
535
0
    return;
536
9.36k
537
9.36k
  partial = sol->partial;
538
9.36k
  if (!partial)
539
2.02k
    return;
540
7.33k
541
7.33k
  if (partial->level == 0 && 
sol->level == 03.11k
) {
542
6.23k
    for (partial = sol->partial; partial; 
partial = sol->partial3.11k
)
543
3.11k
      sol_pop_one(sol);
544
3.11k
    return;
545
3.11k
  }
546
4.22k
547
4.22k
  if (partial->level <= sol->level)
548
8
    return;
549
4.21k
550
4.21k
  if (partial->next && 
partial->next->level == partial->level2.27k
) {
551
2.10k
    if (combine_initial_if_equal(sol) < 0)
552
0
      goto error;
553
2.11k
  } else
554
2.11k
    sol_pop_one(sol);
555
4.21k
556
4.21k
  if (sol->level == 0) {
557
2.22k
    for (partial = sol->partial; partial; 
partial = sol->partial247
)
558
1.98k
      
sol_pop_one(sol)247
;
559
1.98k
    return;
560
1.98k
  }
561
2.23k
562
2.23k
  if (0)
563
0
error:    sol->error = 1;
564
9.36k
}
565
566
static void sol_dec_level(struct isl_sol *sol)
567
2.44k
{
568
2.44k
  if (sol->error)
569
0
    return;
570
2.44k
571
2.44k
  sol->level--;
572
2.44k
573
2.44k
  sol_pop(sol);
574
2.44k
}
575
576
static isl_stat sol_dec_level_wrap(struct isl_tab_callback *cb)
577
2.44k
{
578
2.44k
  struct isl_sol_callback *callback = (struct isl_sol_callback *)cb;
579
2.44k
580
2.44k
  sol_dec_level(callback->sol);
581
2.44k
582
2.44k
  return callback->sol->error ? 
isl_stat_error0
: isl_stat_ok;
583
2.44k
}
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.6k
{
593
21.6k
  struct isl_tab *tab;
594
21.6k
595
21.6k
  if (sol->error)
596
0
    return;
597
21.6k
598
21.6k
  sol->level++;
599
21.6k
  tab = sol->context->op->peek_tab(sol->context);
600
21.6k
  if (isl_tab_push_callback(tab, &sol->dec_level.callback) < 0)
601
0
    sol->error = 1;
602
21.6k
}
603
604
static void scale_rows(struct isl_mat *mat, isl_int m, int n_row)
605
16.3k
{
606
16.3k
  int i;
607
16.3k
608
16.3k
  if (isl_int_is_one(m))
609
16.3k
    
return16.3k
;
610
18
611
51
  
for (i = 0; 18
i < n_row;
++i33
)
612
33
    isl_seq_scale(mat->row[i], mat->row[i], m, mat->n_col);
613
16.3k
}
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
28.5k
{
666
28.5k
  struct isl_basic_set *bset = NULL;
667
28.5k
  struct isl_mat *mat = NULL;
668
28.5k
  unsigned off;
669
28.5k
  int row;
670
28.5k
  isl_int m;
671
28.5k
672
28.5k
  if (sol->error || !tab)
673
0
    goto error;
674
28.5k
675
28.5k
  if (tab->empty && 
!sol->add_empty21.4k
)
676
2.75k
    return;
677
25.8k
  if (sol->context->op->is_empty(sol->context))
678
16.2k
    return;
679
9.58k
680
9.58k
  bset = sol_domain(sol);
681
9.58k
682
9.58k
  if (tab->empty) {
683
2.42k
    sol_push_sol(sol, bset, NULL);
684
2.42k
    return;
685
2.42k
  }
686
7.16k
687
7.16k
  off = 2 + tab->M;
688
7.16k
689
7.16k
  mat = isl_mat_alloc(tab->mat->ctx, 1 + sol->n_out,
690
7.16k
              1 + tab->n_param + tab->n_div);
691
7.16k
  if (!mat)
692
0
    goto error;
693
7.16k
694
7.16k
  isl_int_init(m);
695
7.16k
696
7.16k
  isl_seq_clr(mat->row[0] + 1, mat->n_col - 1);
697
7.16k
  isl_int_set_si(mat->row[0][0], 1);
698
23.4k
  for (row = 0; row < sol->n_out; 
++row16.3k
) {
699
16.3k
    int i = tab->n_param + row;
700
16.3k
    int r, j;
701
16.3k
702
16.3k
    isl_seq_clr(mat->row[1 + row], mat->n_col);
703
16.3k
    if (!tab->var[i].is_row) {
704
0
      if (tab->M)
705
0
        isl_die(mat->ctx, isl_error_invalid,
706
0
          "unbounded optimum", goto error2);
707
0
      continue;
708
16.3k
    }
709
16.3k
710
16.3k
    r = tab->var[i].index;
711
16.3k
    if (tab->M &&
712
16.3k
        isl_int_ne(tab->mat->row[r][2], tab->mat->row[r][0]))
713
16.3k
      
isl_die0
(mat->ctx, isl_error_invalid,
714
16.3k
        "unbounded optimum", goto error2);
715
16.3k
    isl_int_gcd(m, mat->row[0][0], tab->mat->row[r][0]);
716
16.3k
    isl_int_divexact(m, tab->mat->row[r][0], m);
717
16.3k
    scale_rows(mat, m, 1 + row);
718
16.3k
    isl_int_divexact(m, mat->row[0][0], tab->mat->row[r][0]);
719
16.3k
    isl_int_mul(mat->row[1 + row][0], m, tab->mat->row[r][1]);
720
86.8k
    for (j = 0; j < tab->n_param; 
++j70.4k
) {
721
70.4k
      int col;
722
70.4k
      if (tab->var[j].is_row)
723
36.1k
        continue;
724
34.3k
      col = tab->var[j].index;
725
34.3k
      isl_int_mul(mat->row[1 + row][1 + j], m,
726
70.4k
            tab->mat->row[r][off + col]);
727
70.4k
    }
728
18.1k
    for (j = 0; j < tab->n_div; 
++j1.86k
) {
729
1.86k
      int col;
730
1.86k
      if (tab->var[tab->n_var - tab->n_div+j].is_row)
731
173
        continue;
732
1.69k
      col = tab->var[tab->n_var - tab->n_div+j].index;
733
1.69k
      isl_int_mul(mat->row[1 + row][1 + tab->n_param + j], m,
734
1.86k
            tab->mat->row[r][off + col]);
735
1.86k
    }
736
16.3k
    if (sol->max)
737
12.4k
      isl_seq_neg(mat->row[1 + row], mat->row[1 + row],
738
12.4k
            mat->n_col);
739
16.3k
  }
740
7.16k
741
7.16k
  isl_int_clear(m);
742
7.16k
743
7.16k
  sol_push_sol_mat(sol, bset, mat);
744
7.16k
  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.48k
{
761
3.48k
  struct isl_sol_map *sol_map = (struct isl_sol_map *) sol;
762
3.48k
  isl_map_free(sol_map->map);
763
3.48k
  isl_set_free(sol_map->empty);
764
3.48k
}
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.18k
{
773
2.18k
  if (!bset || !sol->empty)
774
0
    goto error;
775
2.18k
776
2.18k
  sol->empty = isl_set_grow(sol->empty, 1);
777
2.18k
  bset = isl_basic_set_simplify(bset);
778
2.18k
  bset = isl_basic_set_finalize(bset);
779
2.18k
  sol->empty = isl_set_add_basic_set(sol->empty, isl_basic_set_copy(bset));
780
2.18k
  if (!sol->empty)
781
0
    goto error;
782
2.18k
  isl_basic_set_free(bset);
783
2.18k
  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.18k
{
792
2.18k
  sol_map_add_empty((struct isl_sol_map *)sol, bset);
793
2.18k
}
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.31k
{
802
3.31k
  isl_basic_map *bmap;
803
3.31k
804
3.31k
  if (sol->sol.error || !dom || !ma)
805
0
    goto error;
806
3.31k
807
3.31k
  bmap = isl_basic_map_from_multi_aff2(ma, sol->sol.rational);
808
3.31k
  bmap = isl_basic_map_intersect_domain(bmap, dom);
809
3.31k
  sol->map = isl_map_grow(sol->map, 1);
810
3.31k
  sol->map = isl_map_add_basic_map(sol->map, bmap);
811
3.31k
  if (!sol->map)
812
0
    sol->sol.error = 1;
813
3.31k
  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.31k
{
823
3.31k
  sol_map_add((struct isl_sol_map *)sol, dom, ma);
824
3.31k
}
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
27.2k
{
836
27.2k
  int i;
837
27.2k
  unsigned off = 2 + tab->M;
838
27.2k
839
27.2k
  isl_int_set(line[0], tab->mat->row[row][1]);
840
158k
  for (i = 0; i < tab->n_param; 
++i131k
) {
841
131k
    if (tab->var[i].is_row)
842
131k
      
isl_int_set_si65.3k
(line[1 + i], 0);
843
131k
    else {
844
66.0k
      int col = tab->var[i].index;
845
66.0k
      isl_int_set(line[1 + i], tab->mat->row[row][off + col]);
846
66.0k
    }
847
131k
  }
848
33.5k
  for (i = 0; i < tab->n_div; 
++i6.38k
) {
849
6.38k
    if (tab->var[tab->n_var - tab->n_div + i].is_row)
850
6.38k
      
isl_int_set_si2.02k
(line[1 + tab->n_param + i], 0);
851
6.38k
    else {
852
4.36k
      int col = tab->var[tab->n_var - tab->n_div + i].index;
853
4.36k
      isl_int_set(line[1 + tab->n_param + i],
854
4.36k
            tab->mat->row[row][off + col]);
855
4.36k
    }
856
6.38k
  }
857
27.2k
}
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
67.1k
{
867
67.1k
  int i;
868
67.1k
  unsigned off = 2 + tab->M;
869
67.1k
870
67.1k
  if (isl_int_ne(tab->mat->row[row1][1], tab->mat->row[row2][1]))
871
67.1k
    
return 061.7k
;
872
5.44k
873
5.44k
  if (tab->M && isl_int_ne(tab->mat->row[row1][2],
874
5.44k
         tab->mat->row[row2][2]))
875
5.44k
    
return 02.19k
;
876
3.24k
877
7.85k
  
for (i = 0; 3.24k
i < tab->n_param + tab->n_div;
++i4.61k
) {
878
7.81k
    int pos = i < tab->n_param ? 
i7.70k
:
879
7.81k
      
tab->n_var - tab->n_div + i - tab->n_param109
;
880
7.81k
    int col;
881
7.81k
882
7.81k
    if (tab->var[pos].is_row)
883
2.88k
      continue;
884
4.92k
    col = tab->var[pos].index;
885
4.92k
    if (isl_int_ne(tab->mat->row[row1][off + col],
886
4.92k
             tab->mat->row[row2][off + col]))
887
4.92k
      
return 03.19k
;
888
7.81k
  }
889
3.24k
  
return 143
;
890
67.1k
}
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
17.2k
{
899
17.2k
  struct isl_vec *ineq;
900
17.2k
901
17.2k
  ineq = isl_vec_alloc(tab->mat->ctx, 1 + tab->n_param + tab->n_div);
902
17.2k
  if (!ineq)
903
0
    return NULL;
904
17.2k
905
17.2k
  get_row_parameter_line(tab, row, ineq->el);
906
17.2k
  if (ineq)
907
17.2k
    ineq = isl_vec_normalize(ineq);
908
17.2k
909
17.2k
  return ineq;
910
17.2k
}
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
607
{
922
607
  isl_ctx *ctx = isl_vec_get_ctx(div);
923
607
  int len = div->size - 2;
924
607
925
607
  isl_seq_gcd(div->el + 2, len, &ctx->normalize_gcd);
926
607
  isl_int_gcd(ctx->normalize_gcd, ctx->normalize_gcd, div->el[0]);
927
607
928
607
  if (isl_int_is_one(ctx->normalize_gcd))
929
607
    
return558
;
930
49
931
49
  isl_int_divexact(div->el[0], div->el[0], ctx->normalize_gcd);
932
49
  isl_int_fdiv_q(div->el[1], div->el[1], ctx->normalize_gcd);
933
607
  isl_seq_scale_down(div->el + 2, div->el + 2, ctx->normalize_gcd, len);
934
607
}
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
517
{
949
517
  struct isl_vec *div;
950
517
951
517
  div = isl_vec_alloc(tab->mat->ctx, 1 + 1 + tab->n_param + tab->n_div);
952
517
  if (!div)
953
0
    return NULL;
954
517
955
517
  isl_int_set(div->el[0], tab->mat->row[row][0]);
956
517
  get_row_parameter_line(tab, row, div->el + 1);
957
517
  isl_seq_neg(div->el + 1, div->el + 1, div->size - 1);
958
517
  normalize_div(div);
959
517
  isl_seq_fdiv_r(div->el + 1, div->el + 1, div->el[0], div->size - 1);
960
517
961
517
  return div;
962
517
}
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
90
{
977
90
  struct isl_vec *div;
978
90
979
90
  div = isl_vec_alloc(tab->mat->ctx, 1 + 1 + tab->n_param + tab->n_div);
980
90
  if (!div)
981
0
    return NULL;
982
90
983
90
  isl_int_set(div->el[0], tab->mat->row[row][0]);
984
90
  get_row_parameter_line(tab, row, div->el + 1);
985
90
  normalize_div(div);
986
90
  isl_seq_fdiv_r(div->el + 1, div->el + 1, div->el[0], div->size - 1);
987
90
988
90
  return div;
989
90
}
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
90
{
1004
90
  unsigned total;
1005
90
  unsigned div_pos;
1006
90
  struct isl_vec *ineq;
1007
90
1008
90
  if (!bset)
1009
0
    return NULL;
1010
90
1011
90
  total = isl_basic_set_total_dim(bset);
1012
90
  div_pos = 1 + total - bset->n_div + div;
1013
90
1014
90
  ineq = isl_vec_alloc(bset->ctx, 1 + total);
1015
90
  if (!ineq)
1016
0
    return NULL;
1017
90
1018
90
  isl_seq_cpy(ineq->el, bset->div[div] + 1, 1 + total);
1019
90
  isl_int_neg(ineq->el[div_pos], bset->div[div][0]);
1020
90
  return ineq;
1021
90
}
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
90
{
1038
90
  isl_seq_fdiv_q(tab->mat->row[row] + 1, tab->mat->row[row] + 1,
1039
90
      tab->mat->row[row][0], 1 + tab->M + tab->n_col);
1040
90
1041
90
  isl_int_set_si(tab->mat->row[row][0], 1);
1042
90
1043
90
  if (tab->var[tab->n_var - tab->n_div + div].is_row) {
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
90
  } else {
1053
90
    int dcol = tab->var[tab->n_var - tab->n_div + div].index;
1054
90
1055
90
    isl_int_add_ui(tab->mat->row[row][2 + tab->M + dcol],
1056
90
        tab->mat->row[row][2 + tab->M + dcol], 1);
1057
90
  }
1058
90
1059
90
  return tab;
1060
0
error:
1061
0
  isl_tab_free(tab);
1062
0
  return NULL;
1063
90
}
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
199k
{
1075
199k
  int i;
1076
199k
  int col;
1077
199k
  unsigned off = 2 + tab->M;
1078
199k
1079
199k
  if (tab->M) {
1080
101k
    if (isl_int_is_pos(tab->mat->row[row][2]))
1081
101k
      
return 053.2k
;
1082
47.9k
    if (isl_int_is_neg(tab->mat->row[row][2]))
1083
47.9k
      
return 10
;
1084
146k
  }
1085
146k
1086
146k
  if (isl_int_is_nonneg(tab->mat->row[row][1]))
1087
146k
    
return 0132k
;
1088
28.0k
  
for (i = 0; 14.2k
i < tab->n_param;
++i13.8k
) {
1089
24.1k
    /* Eliminated parameter */
1090
24.1k
    if (tab->var[i].is_row)
1091
7.38k
      continue;
1092
16.7k
    col = tab->var[i].index;
1093
16.7k
    if (isl_int_is_zero(tab->mat->row[row][off + col]))
1094
16.7k
      
continue6.41k
;
1095
10.3k
    if (!tab->var[i].is_nonneg)
1096
9.36k
      return 0;
1097
960
    if (isl_int_is_pos(tab->mat->row[row][off + col]))
1098
960
      
return 0919
;
1099
24.1k
  }
1100
14.2k
  
for (i = 0; 3.96k
i < tab->n_div4.42k
;
++i452
) {
1101
630
    if (tab->var[tab->n_var - tab->n_div + i].is_row)
1102
404
      continue;
1103
226
    col = tab->var[tab->n_var - tab->n_div + i].index;
1104
226
    if (isl_int_is_zero(tab->mat->row[row][off + col]))
1105
226
      
continue48
;
1106
178
    if (!tab->var[tab->n_var - tab->n_div + i].is_nonneg)
1107
151
      return 0;
1108
27
    if (isl_int_is_pos(tab->mat->row[row][off + col]))
1109
27
      return 0;
1110
630
  }
1111
3.96k
  
return 13.79k
;
1112
199k
}
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
32.4k
{
1124
32.4k
  int i;
1125
32.4k
  int col;
1126
32.4k
  unsigned off = 2 + tab->M;
1127
32.4k
1128
32.4k
  if (tab->M) {
1129
32.4k
    if (isl_int_is_pos(tab->mat->row[row][2]))
1130
32.4k
      
return 113.6k
;
1131
18.7k
    if (isl_int_is_neg(tab->mat->row[row][2]))
1132
18.7k
      
return 00
;
1133
18.7k
  }
1134
18.7k
1135
18.7k
  if (isl_int_is_neg(tab->mat->row[row][1]))
1136
18.7k
    
return 06.08k
;
1137
48.2k
  
for (i = 0; 12.7k
i < tab->n_param;
++i35.5k
) {
1138
41.4k
    /* Eliminated parameter */
1139
41.4k
    if (tab->var[i].is_row)
1140
17.2k
      continue;
1141
24.2k
    col = tab->var[i].index;
1142
24.2k
    if (isl_int_is_zero(tab->mat->row[row][off + col]))
1143
24.2k
      
continue15.2k
;
1144
8.95k
    if (!tab->var[i].is_nonneg)
1145
1.84k
      return 0;
1146
7.10k
    if (isl_int_is_neg(tab->mat->row[row][off + col]))
1147
7.10k
      
return 04.04k
;
1148
41.4k
  }
1149
12.7k
  
for (i = 0; 6.81k
i < tab->n_div11.1k
;
++i4.38k
) {
1150
4.56k
    if (tab->var[tab->n_var - tab->n_div + i].is_row)
1151
3.75k
      continue;
1152
810
    col = tab->var[tab->n_var - tab->n_div + i].index;
1153
810
    if (isl_int_is_zero(tab->mat->row[row][off + col]))
1154
810
      
continue616
;
1155
194
    if (!tab->var[tab->n_var - tab->n_div + i].is_nonneg)
1156
179
      return 0;
1157
15
    if (isl_int_is_neg(tab->mat->row[row][off + col]))
1158
15
      
return 07
;
1159
4.56k
  }
1160
6.81k
  
return 16.62k
;
1161
32.4k
}
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.96k
{
1182
4.96k
  int i;
1183
4.96k
  isl_int *tr;
1184
4.96k
1185
4.96k
  tr = tab->mat->row[row] + 2 + tab->M;
1186
4.96k
1187
21.1k
  for (i = tab->n_param; i < tab->n_var - tab->n_div; 
++i16.1k
) {
1188
21.1k
    int s1, s2;
1189
21.1k
    isl_int *r;
1190
21.1k
1191
21.1k
    if (!tab->var[i].is_row) {
1192
7.20k
      if (tab->var[i].index == col1)
1193
701
        return col2;
1194
6.50k
      if (tab->var[i].index == col2)
1195
464
        return col1;
1196
6.03k
      continue;
1197
6.03k
    }
1198
13.9k
1199
13.9k
    if (tab->var[i].index == row)
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 && 
s2 == 010.4k
)
1206
8.96k
      continue;
1207
4.83k
    if (s1 < s2)
1208
1.67k
      return col1;
1209
3.16k
    if (s2 < s1)
1210
887
      return col2;
1211
2.27k
1212
2.27k
    isl_int_mul(tmp, r[col2], tr[col1]);
1213
2.27k
    isl_int_submul(tmp, r[col1], tr[col2]);
1214
2.27k
    if (isl_int_is_pos(tmp))
1215
2.27k
      
return col1836
;
1216
1.43k
    if (isl_int_is_neg(tmp))
1217
1.43k
      
return col2410
;
1218
21.1k
  }
1219
4.96k
  
return -10
;
1220
4.96k
}
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
12.2k
{
1230
12.2k
  int j;
1231
12.2k
  int col = tab->n_col;
1232
12.2k
  isl_int *tr;
1233
12.2k
  isl_int tmp;
1234
12.2k
1235
12.2k
  tr = tab->mat->row[row] + 2 + tab->M;
1236
12.2k
1237
12.2k
  isl_int_init(tmp);
1238
12.2k
1239
90.6k
  for (j = tab->n_dead; j < tab->n_col; 
++j78.4k
) {
1240
78.4k
    if (tab->col_var[j] >= 0 &&
1241
78.4k
        
(59.2k
tab->col_var[j] < tab->n_param59.2k
||
1242
59.2k
        
tab->col_var[j] >= tab->n_var - tab->n_div41.3k
))
1243
19.8k
      continue;
1244
58.5k
1245
58.5k
    if (!isl_int_is_pos(tr[j]))
1246
58.5k
      
continue44.5k
;
1247
14.0k
1248
14.0k
    if (col == tab->n_col)
1249
9.03k
      col = j;
1250
4.96k
    else
1251
4.96k
      col = lexmin_col_pair(tab, row, col, j, tmp);
1252
14.0k
    isl_assert(tab->mat->ctx, col >= 0, goto error);
1253
78.4k
  }
1254
12.2k
1255
12.2k
  isl_int_clear(tmp);
1256
12.2k
  return col;
1257
0
error:
1258
0
  isl_int_clear(tmp);
1259
0
  return -1;
1260
12.2k
}
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
42.9k
{
1271
42.9k
  int row;
1272
42.9k
1273
42.9k
  if (tab->M)
1274
261k
    
for (row = tab->n_redundant; 34.6k
row < tab->n_row;
++row227k
) {
1275
231k
      if (!isl_tab_var_from_row(tab, row)->is_nonneg)
1276
55.0k
        continue;
1277
176k
      if (!isl_int_is_neg(tab->mat->row[row][2]))
1278
176k
        
continue172k
;
1279
4.49k
      if (tab->row_sign)
1280
4.49k
        tab->row_sign[row] = isl_tab_row_neg;
1281
34.6k
      return row;
1282
34.6k
    }
1283
329k
  
for (row = tab->n_redundant; 38.4k
row < tab->n_row;
++row290k
) {
1284
298k
    if (!isl_tab_var_from_row(tab, row)->is_nonneg)
1285
47.4k
      continue;
1286
251k
    if (tab->row_sign) {
1287
152k
      if (tab->row_sign[row] == 0 &&
1288
152k
          
is_obviously_neg(tab, row)101k
)
1289
270
        tab->row_sign[row] = isl_tab_row_neg;
1290
152k
      if (tab->row_sign[row] != isl_tab_row_neg)
1291
148k
        continue;
1292
98.4k
    } else if (!is_obviously_neg(tab, row))
1293
94.8k
      continue;
1294
7.71k
    return row;
1295
7.71k
  }
1296
38.4k
  
return -130.7k
;
1297
42.9k
}
1298
1299
/* Check whether the invariant that all columns are lexico-positive
1300
 * is satisfied.  This function is not called from the current code
1301
 * but is useful during debugging.
1302
 */
1303
static void check_lexpos(struct isl_tab *tab) __attribute__ ((unused));
1304
static void check_lexpos(struct isl_tab *tab)
1305
0
{
1306
0
  unsigned off = 2 + tab->M;
1307
0
  int col;
1308
0
  int var;
1309
0
  int row;
1310
0
1311
0
  for (col = tab->n_dead; col < tab->n_col; ++col) {
1312
0
    if (tab->col_var[col] >= 0 &&
1313
0
        (tab->col_var[col] < tab->n_param ||
1314
0
         tab->col_var[col] >= tab->n_var - tab->n_div))
1315
0
      continue;
1316
0
    for (var = tab->n_param; var < tab->n_var - tab->n_div; ++var) {
1317
0
      if (!tab->var[var].is_row) {
1318
0
        if (tab->var[var].index == col)
1319
0
          break;
1320
0
        else
1321
0
          continue;
1322
0
      }
1323
0
      row = tab->var[var].index;
1324
0
      if (isl_int_is_zero(tab->mat->row[row][off + col]))
1325
0
        continue;
1326
0
      if (isl_int_is_pos(tab->mat->row[row][off + col]))
1327
0
        break;
1328
0
      fprintf(stderr, "lexneg column %d (row %d)\n",
1329
0
        col, row);
1330
0
    }
1331
0
    if (var >= tab->n_var - tab->n_div)
1332
0
      fprintf(stderr, "zero column %d\n", col);
1333
0
  }
1334
0
}
1335
1336
/* Report to the caller that the given constraint is part of an encountered
1337
 * conflict.
1338
 */
1339
static int report_conflicting_constraint(struct isl_tab *tab, int con)
1340
1.17k
{
1341
1.17k
  return tab->conflict(con, tab->conflict_user);
1342
1.17k
}
1343
1344
/* Given a conflicting row in the tableau, report all constraints
1345
 * involved in the row to the caller.  That is, the row itself
1346
 * (if it represents a constraint) and all constraint columns with
1347
 * non-zero (and therefore negative) coefficients.
1348
 */
1349
static int report_conflict(struct isl_tab *tab, int row)
1350
3.17k
{
1351
3.17k
  int j;
1352
3.17k
  isl_int *tr;
1353
3.17k
1354
3.17k
  if (!tab->conflict)
1355
2.71k
    return 0;
1356
459
1357
459
  if (tab->row_var[row] < 0 &&
1358
459
      report_conflicting_constraint(tab, ~tab->row_var[row]) < 0)
1359
0
    return -1;
1360
459
1361
459
  tr = tab->mat->row[row] + 2 + tab->M;
1362
459
1363
5.37k
  for (j = tab->n_dead; j < tab->n_col; 
++j4.92k
) {
1364
4.92k
    if (tab->col_var[j] >= 0 &&
1365
4.92k
        
(3.99k
tab->col_var[j] < tab->n_param3.99k
||
1366
3.99k
        tab->col_var[j] >= tab->n_var - tab->n_div))
1367
0
      continue;
1368
4.92k
1369
4.92k
    if (!isl_int_is_neg(tr[j]))
1370
4.92k
      
continue4.17k
;
1371
749
1372
749
    if (tab->col_var[j] < 0 &&
1373
749
        
report_conflicting_constraint(tab, ~tab->col_var[j]) < 0713
)
1374
0
      return -1;
1375
4.92k
  }
1376
459
1377
459
  return 0;
1378
3.17k
}
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
33.9k
{
1389
33.9k
  int row, col;
1390
33.9k
1391
33.9k
  if (!tab)
1392
0
    return -1;
1393
33.9k
  if (tab->empty)
1394
26
    return 0;
1395
42.9k
  
while (33.8k
(row = first_neg(tab)) != -1) {
1396
12.2k
    col = lexmin_pivot_col(tab, row);
1397
12.2k
    if (col >= tab->n_col) {
1398
3.17k
      if (report_conflict(tab, row) < 0)
1399
0
        return -1;
1400
3.17k
      if (isl_tab_mark_empty(tab) < 0)
1401
0
        return -1;
1402
3.17k
      return 0;
1403
3.17k
    }
1404
9.03k
    if (col < 0)
1405
0
      return -1;
1406
9.03k
    if (isl_tab_pivot(tab, row, col) < 0)
1407
0
      return -1;
1408
12.2k
  }
1409
33.8k
  
return 030.7k
;
1410
33.9k
}
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
23.0k
{
1435
23.0k
  unsigned off = 2 + tab->M;
1436
23.0k
  int i;
1437
23.0k
1438
81.5k
  for (i = tab->n_var - tab->n_div - 1; i >= 0 && 
i >= tab->n_param80.5k
;
--i58.4k
) {
1439
71.2k
    int col;
1440
71.2k
    if (tab->var[i].is_row)
1441
41.8k
      continue;
1442
29.3k
    col = tab->var[i].index;
1443
29.3k
    if (col <= tab->n_dead)
1444
2.26k
      continue;
1445
27.1k
    if (!isl_int_is_zero(tab->mat->row[row][off + col]))
1446
27.1k
      
return col12.7k
;
1447
71.2k
  }
1448
27.2k
  
for (i = tab->n_dead; 10.2k
i < tab->n_col;
++i16.9k
) {
1449
27.1k
    if (isl_int_is_one(tab->mat->row[row][off + i]))
1450
27.1k
      
return i4.78k
;
1451
22.3k
    if (isl_int_is_negone(tab->mat->row[row][off + i]))
1452
22.3k
      
return i5.42k
;
1453
27.1k
  }
1454
10.2k
  
return -183
;
1455
23.0k
}
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
23.0k
{
1468
23.0k
  int i;
1469
23.0k
  int r;
1470
23.0k
1471
23.0k
  if (!tab)
1472
0
    return NULL;
1473
23.0k
  r = isl_tab_add_row(tab, eq);
1474
23.0k
  if (r < 0)
1475
0
    goto error;
1476
23.0k
1477
23.0k
  r = tab->con[r].index;
1478
23.0k
  i = last_var_col_or_int_par_col(tab, r);
1479
23.0k
  if (i < 0) {
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
22.9k
  } else {
1491
22.9k
    if (isl_tab_pivot(tab, r, i) < 0)
1492
0
      goto error;
1493
22.9k
    if (isl_tab_kill_col(tab, i) < 0)
1494
0
      goto error;
1495
22.9k
    tab->n_eq++;
1496
22.9k
  }
1497
23.0k
1498
23.0k
  return tab;
1499
0
error:
1500
0
  isl_tab_free(tab);
1501
0
  return NULL;
1502
23.0k
}
1503
1504
/* Check if the given row is a pure constant.
1505
 */
1506
static int is_constant(struct isl_tab *tab, int row)
1507
305
{
1508
305
  unsigned off = 2 + tab->M;
1509
305
1510
305
  return isl_seq_first_non_zero(tab->mat->row[row] + off + tab->n_dead,
1511
305
          tab->n_col - tab->n_dead) == -1;
1512
305
}
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
305
{
1529
305
  int r1, r2;
1530
305
  int row;
1531
305
  struct isl_tab_undo *snap;
1532
305
1533
305
  if (!tab)
1534
0
    return -1;
1535
305
  snap = isl_tab_snap(tab);
1536
305
  r1 = isl_tab_add_row(tab, eq);
1537
305
  if (r1 < 0)
1538
0
    return -1;
1539
305
  tab->con[r1].is_nonneg = 1;
1540
305
  if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r1]) < 0)
1541
0
    return -1;
1542
305
1543
305
  row = tab->con[r1].index;
1544
305
  if (is_constant(tab, row)) {
1545
54
    if (!isl_int_is_zero(tab->mat->row[row][1]) ||
1546
54
        (tab->M && 
!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 (isl_tab_rollback(tab, snap) < 0)
1552
0
      return -1;
1553
54
    return 0;
1554
54
  }
1555
251
1556
251
  if (restore_lexmin(tab) < 0)
1557
0
    return -1;
1558
251
  if (tab->empty)
1559
14
    return 0;
1560
237
1561
237
  isl_seq_neg(eq, eq, 1 + tab->n_var);
1562
237
1563
237
  r2 = isl_tab_add_row(tab, eq);
1564
237
  if (r2 < 0)
1565
0
    return -1;
1566
237
  tab->con[r2].is_nonneg = 1;
1567
237
  if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r2]) < 0)
1568
0
    return -1;
1569
237
1570
237
  if (restore_lexmin(tab) < 0)
1571
0
    return -1;
1572
237
  if (tab->empty)
1573
0
    return 0;
1574
237
1575
237
  if (!tab->con[r1].is_row) {
1576
0
    if (isl_tab_kill_col(tab, tab->con[r1].index) < 0)
1577
0
      return -1;
1578
237
  } else if (!tab->con[r2].is_row) {
1579
0
    if (isl_tab_kill_col(tab, tab->con[r2].index) < 0)
1580
0
      return -1;
1581
237
  }
1582
237
1583
237
  if (tab->bmap) {
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 (!tab->bmap)
1593
0
      return -1;
1594
237
  }
1595
237
1596
237
  return 0;
1597
237
}
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
22.9k
{
1607
22.9k
  int r;
1608
22.9k
1609
22.9k
  if (!tab)
1610
0
    return NULL;
1611
22.9k
  if (tab->bmap) {
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 (!tab->bmap)
1616
0
      goto error;
1617
22.9k
  }
1618
22.9k
  r = isl_tab_add_row(tab, ineq);
1619
22.9k
  if (r < 0)
1620
0
    goto error;
1621
22.9k
  tab->con[r].is_nonneg = 1;
1622
22.9k
  if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r]) < 0)
1623
0
    goto error;
1624
22.9k
  if (isl_tab_row_is_redundant(tab, tab->con[r].index)) {
1625
41
    if (isl_tab_mark_redundant(tab, tab->con[r].index) < 0)
1626
0
      goto error;
1627
41
    return tab;
1628
41
  }
1629
22.8k
1630
22.8k
  if (restore_lexmin(tab) < 0)
1631
0
    goto error;
1632
22.8k
  if (!tab->empty && 
tab->con[r].is_row22.4k
&&
1633
22.8k
     
isl_tab_row_is_redundant(tab, tab->con[r].index)17.5k
)
1634
0
    if (isl_tab_mark_redundant(tab, tab->con[r].index) < 0)
1635
0
      goto error;
1636
22.8k
  return tab;
1637
0
error:
1638
0
  isl_tab_free(tab);
1639
0
  return NULL;
1640
22.9k
}
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
22.5k
{
1646
22.5k
  int i;
1647
22.5k
  int col;
1648
22.5k
  unsigned off = 2 + tab->M;
1649
22.5k
1650
100k
  for (i = 0; i < tab->n_param; 
++i77.6k
) {
1651
78.2k
    /* Eliminated parameter */
1652
78.2k
    if (tab->var[i].is_row)
1653
37.3k
      continue;
1654
40.8k
    col = tab->var[i].index;
1655
40.8k
    if (!isl_int_is_divisible_by(tab->mat->row[row][off + col],
1656
40.8k
            tab->mat->row[row][0]))
1657
40.8k
      
return 0531
;
1658
78.2k
  }
1659
27.6k
  
for (i = 0; 22.0k
i < tab->n_div;
++i5.54k
) {
1660
5.62k
    if (tab->var[tab->n_var - tab->n_div + i].is_row)
1661
1.38k
      continue;
1662
4.23k
    col = tab->var[tab->n_var - tab->n_div + i].index;
1663
4.23k
    if (!isl_int_is_divisible_by(tab->mat->row[row][off + col],
1664
4.23k
            tab->mat->row[row][0]))
1665
4.23k
      
return 076
;
1666
5.62k
  }
1667
22.0k
  
return 121.9k
;
1668
22.5k
}
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.09k
{
1674
1.09k
  int i;
1675
1.09k
  unsigned off = 2 + tab->M;
1676
1.09k
1677
3.07k
  for (i = tab->n_dead; i < tab->n_col; 
++i1.97k
) {
1678
2.98k
    if (tab->col_var[i] >= 0 &&
1679
2.98k
        
(1.78k
tab->col_var[i] < tab->n_param1.78k
||
1680
1.78k
         
tab->col_var[i] >= tab->n_var - tab->n_div142
))
1681
1.75k
      continue;
1682
1.22k
    if (!isl_int_is_divisible_by(tab->mat->row[row][off + i],
1683
1.22k
            tab->mat->row[row][0]))
1684
1.22k
      
return 01.00k
;
1685
2.98k
  }
1686
1.09k
  
return 190
;
1687
1.09k
}
1688
1689
/* Check if the constant term is integral.
1690
 */
1691
static int integer_constant(struct isl_tab *tab, int row)
1692
22.5k
{
1693
22.5k
  return isl_int_is_divisible_by(tab->mat->row[row][1],
1694
22.5k
          tab->mat->row[row][0]);
1695
22.5k
}
1696
1697
#define I_CST 1 << 0
1698
#define I_PAR 1 << 1
1699
#define I_VAR 1 << 2
1700
1701
/* Check for next (non-parameter) variable after "var" (first if var == -1)
1702
 * that is non-integer and therefore requires a cut and return
1703
 * the index of the variable.
1704
 * For parametric tableaus, there are three parts in a row,
1705
 * the constant, the coefficients of the parameters and the rest.
1706
 * For each part, we check whether the coefficients in that part
1707
 * are all integral and if so, set the corresponding flag in *f.
1708
 * If the constant and the parameter part are integral, then the
1709
 * current sample value is integral and no cut is required
1710
 * (irrespective of whether the variable part is integral).
1711
 */
1712
static int next_non_integer_var(struct isl_tab *tab, int var, int *f)
1713
8.98k
{
1714
8.98k
  var = var < 0 ? tab->n_param : 
var + 10
;
1715
8.98k
1716
36.9k
  for (; var < tab->n_var - tab->n_div; 
++var27.9k
) {
1717
29.0k
    int flags = 0;
1718
29.0k
    int row;
1719
29.0k
    if (!tab->var[var].is_row)
1720
6.43k
      continue;
1721
22.5k
    row = tab->var[var].index;
1722
22.5k
    if (integer_constant(tab, row))
1723
22.5k
      
ISL_FL_SET21.7k
(flags, I_CST);
1724
22.5k
    if (integer_parameter(tab, row))
1725
22.5k
      
ISL_FL_SET21.9k
(flags, I_PAR);
1726
22.5k
    if (ISL_FL_ISSET(flags, I_CST) && 
ISL_FL_ISSET21.7k
(flags, I_PAR))
1727
22.5k
      
continue21.5k
;
1728
1.09k
    if (integer_variable(tab, row))
1729
1.09k
      
ISL_FL_SET90
(flags, I_VAR);
1730
29.0k
    *f = flags;
1731
29.0k
    return var;
1732
29.0k
  }
1733
8.98k
  
return -17.88k
;
1734
8.98k
}
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
8.20k
{
1748
8.20k
  int var = next_non_integer_var(tab, -1, f);
1749
8.20k
1750
8.20k
  return var < 0 ? 
-17.16k
:
tab->var[var].index1.04k
;
1751
8.20k
}
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
488
{
1777
488
  int i;
1778
488
  int r;
1779
488
  isl_int *r_row;
1780
488
  unsigned off = 2 + tab->M;
1781
488
1782
488
  if (isl_tab_extend_cons(tab, 1) < 0)
1783
0
    return -1;
1784
488
  r = isl_tab_allocate_con(tab);
1785
488
  if (r < 0)
1786
0
    return -1;
1787
488
1788
488
  r_row = tab->mat->row[tab->con[r].index];
1789
488
  isl_int_set(r_row[0], tab->mat->row[row][0]);
1790
488
  isl_int_neg(r_row[1], tab->mat->row[row][1]);
1791
488
  isl_int_fdiv_r(r_row[1], r_row[1], tab->mat->row[row][0]);
1792
488
  isl_int_neg(r_row[1], r_row[1]);
1793
488
  if (tab->M)
1794
488
    
isl_int_set_si433
(r_row[2], 0);
1795
4.55k
  for (i = 0; i < tab->n_col; 
++i4.06k
)
1796
4.06k
    isl_int_fdiv_r(r_row[off + i],
1797
488
      tab->mat->row[row][off + i], tab->mat->row[row][0]);
1798
488
1799
488
  tab->con[r].is_nonneg = 1;
1800
488
  if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r]) < 0)
1801
0
    return -1;
1802
488
  if (tab->row_sign)
1803
433
    tab->row_sign[tab->con[r].index] = isl_tab_row_neg;
1804
488
1805
488
  return tab->con[r].index;
1806
488
}
1807
1808
0
#define CUT_ALL 1
1809
1.25k
#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.20k
{
1829
1.20k
  int var;
1830
1.20k
  int row;
1831
1.20k
  int flags;
1832
1.20k
1833
1.20k
  if (!tab)
1834
0
    return NULL;
1835
1.20k
  if (tab->empty)
1836
469
    return tab;
1837
731
1838
784
  
while (731
(var = next_non_integer_var(tab, -1, &flags)) != -1) {
1839
55
    do {
1840
55
      if (ISL_FL_ISSET(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 (cutting_strategy == CUT_ONE)
1850
55
        break;
1851
0
    } while ((var = next_non_integer_var(tab, var, &flags)) != -1);
1852
55
    if (restore_lexmin(tab) < 0)
1853
0
      goto error;
1854
55
    if (tab->empty)
1855
2
      break;
1856
55
  }
1857
731
  return tab;
1858
0
error:
1859
0
  isl_tab_free(tab);
1860
0
  return NULL;
1861
1.20k
}
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
15.3k
{
1869
15.3k
  int i;
1870
15.3k
  isl_int v;
1871
15.3k
1872
15.3k
  if (!tab)
1873
0
    return NULL;
1874
15.3k
1875
15.3k
  isl_assert(tab->mat->ctx, tab->bmap, goto error);
1876
15.3k
  isl_assert(tab->mat->ctx, tab->samples, goto error);
1877
15.3k
  isl_assert(tab->mat->ctx, tab->samples->n_col == 1 + tab->n_var, goto error);
1878
15.3k
1879
15.3k
  isl_int_init(v);
1880
38.1k
  for (i = tab->n_outside; i < tab->n_sample; 
++i22.8k
) {
1881
22.8k
    int sgn;
1882
22.8k
    isl_seq_inner_product(ineq, tab->samples->row[i],
1883
22.8k
          1 + tab->n_var, &v);
1884
22.8k
    sgn = isl_int_sgn(v);
1885
22.8k
    if (eq ? 
(sgn == 0)9.41k
:
(sgn >= 0)13.4k
)
1886
16.0k
      continue;
1887
6.77k
    tab = isl_tab_drop_sample(tab, i);
1888
6.77k
    if (!tab)
1889
0
      break;
1890
22.8k
  }
1891
15.3k
  isl_int_clear(v);
1892
15.3k
1893
15.3k
  return tab;
1894
0
error:
1895
0
  isl_tab_free(tab);
1896
0
  return NULL;
1897
15.3k
}
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; i < tab->n_var; ++i) {
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, CUT_ALL);
1942
0
  if (!tab)
1943
0
    goto error;
1944
0
1945
0
  if (!tab->empty && sample_is_finite(tab)) {
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 (!tab->empty && isl_tab_rollback(tab, snap) < 0)
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
28.0k
{
1968
28.0k
  int i;
1969
28.0k
  isl_int v;
1970
28.0k
1971
28.0k
  if (!tab)
1972
0
    return -1;
1973
28.0k
1974
28.0k
  isl_assert(tab->mat->ctx, tab->bmap, return -1);
1975
28.0k
  isl_assert(tab->mat->ctx, tab->samples, return -1);
1976
28.0k
  isl_assert(tab->mat->ctx, tab->samples->n_col == 1 + tab->n_var, return -1);
1977
28.0k
1978
28.0k
  isl_int_init(v);
1979
46.7k
  for (i = tab->n_outside; i < tab->n_sample; 
++i18.6k
) {
1980
28.0k
    int sgn;
1981
28.0k
    isl_seq_inner_product(ineq, tab->samples->row[i],
1982
28.0k
          1 + tab->n_var, &v);
1983
28.0k
    sgn = isl_int_sgn(v);
1984
28.0k
    if (eq ? 
(sgn == 0)9.30k
:
(sgn >= 0)18.7k
)
1985
9.39k
      break;
1986
28.0k
  }
1987
28.0k
  isl_int_clear(v);
1988
28.0k
1989
28.0k
  return i < tab->n_sample;
1990
28.0k
}
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
580
{
1999
580
  int i;
2000
580
  int r;
2001
580
  struct isl_mat *samples;
2002
580
  int nonneg;
2003
580
2004
580
  r = isl_tab_insert_div(tab, pos, div, add_ineq, user);
2005
580
  if (r < 0)
2006
0
    return isl_bool_error;
2007
580
  nonneg = tab->var[r].is_nonneg;
2008
580
  tab->var[r].frozen = 1;
2009
580
2010
580
  samples = isl_mat_extend(tab->samples,
2011
580
      tab->n_sample, 1 + tab->n_var);
2012
580
  tab->samples = samples;
2013
580
  if (!samples)
2014
0
    return isl_bool_error;
2015
1.55k
  
for (i = tab->n_outside; 580
i < samples->n_row;
++i976
) {
2016
976
    isl_seq_inner_product(div->el + 1, samples->row[i],
2017
976
      div->size - 1, &samples->row[i][samples->n_col - 1]);
2018
976
    isl_int_fdiv_q(samples->row[i][samples->n_col - 1],
2019
976
             samples->row[i][samples->n_col - 1], div->el[0]);
2020
976
  }
2021
580
  tab->samples = isl_mat_move_cols(tab->samples, 1 + pos,
2022
580
          1 + tab->n_var - 1, 1);
2023
580
  if (!tab->samples)
2024
0
    return isl_bool_error;
2025
580
2026
580
  return nonneg;
2027
580
}
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
580
{
2043
580
  int r;
2044
580
  int pos;
2045
580
  isl_bool nonneg;
2046
580
  struct isl_tab *context_tab = context->op->peek_tab(context);
2047
580
2048
580
  if (!tab || !context_tab)
2049
0
    goto error;
2050
580
2051
580
  pos = context_tab->n_var - context->n_unknown;
2052
580
  if ((nonneg = context->op->insert_div(context, pos, div)) < 0)
2053
0
    goto error;
2054
580
2055
580
  if (!context->op->is_ok(context))
2056
0
    goto error;
2057
580
2058
580
  pos = tab->n_var - context->n_unknown;
2059
580
  if (isl_tab_extend_vars(tab, 1) < 0)
2060
0
    goto error;
2061
580
  r = isl_tab_insert_var(tab, pos);
2062
580
  if (r < 0)
2063
0
    goto error;
2064
580
  if (nonneg)
2065
80
    tab->var[r].is_nonneg = 1;
2066
580
  tab->var[r].frozen = 1;
2067
580
  tab->n_div++;
2068
580
2069
580
  return tab->n_div - 1 - context->n_unknown;
2070
0
error:
2071
0
  context->op->invalidate(context);
2072
0
  return -1;
2073
580
}
2074
2075
static int find_div(struct isl_tab *tab, isl_int *div, isl_int denom)
2076
607
{
2077
607
  int i;
2078
607
  unsigned total = isl_basic_map_total_dim(tab->bmap);
2079
607
2080
1.32k
  for (i = 0; i < tab->bmap->n_div; 
++i720
) {
2081
747
    if (isl_int_ne(tab->bmap->div[i][0], denom))
2082
747
      
continue664
;
2083
83
    if (!isl_seq_eq(tab->bmap->div[i] + 1, div, 1 + total))
2084
56
      continue;
2085
27
    return i;
2086
27
  }
2087
607
  
return -1580
;
2088
607
}
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
607
{
2096
607
  int d;
2097
607
  struct isl_tab *context_tab = context->op->peek_tab(context);
2098
607
2099
607
  if (!context_tab)
2100
0
    return -1;
2101
607
2102
607
  d = find_div(context_tab, div->el + 1, div->el[0]);
2103
607
  if (d != -1)
2104
27
    return d;
2105
580
2106
580
  return add_div(tab, context, div);
2107
580
}
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
517
{
2135
517
  struct isl_vec *div;
2136
517
  int d;
2137
517
  int i;
2138
517
  int r;
2139
517
  isl_int *r_row;
2140
517
  int col;
2141
517
  int n;
2142
517
  unsigned off = 2 + tab->M;
2143
517
2144
517
  if (!context)
2145
0
    return -1;
2146
517
2147
517
  div = get_row_parameter_div(tab, row);
2148
517
  if (!div)
2149
0
    return -1;
2150
517
2151
517
  n = tab->n_div - context->n_unknown;
2152
517
  d = context->op->get_div(context, tab, div);
2153
517
  isl_vec_free(div);
2154
517
  if (d < 0)
2155
0
    return -1;
2156
517
2157
517
  if (isl_tab_extend_cons(tab, 1) < 0)
2158
0
    return -1;
2159
517
  r = isl_tab_allocate_con(tab);
2160
517
  if (r < 0)
2161
0
    return -1;
2162
517
2163
517
  r_row = tab->mat->row[tab->con[r].index];
2164
517
  isl_int_set(r_row[0], tab->mat->row[row][0]);
2165
517
  isl_int_neg(r_row[1], tab->mat->row[row][1]);
2166
517
  isl_int_fdiv_r(r_row[1], r_row[1], tab->mat->row[row][0]);
2167
517
  isl_int_neg(r_row[1], r_row[1]);
2168
517
  if (tab->M)
2169
517
    isl_int_set_si(r_row[2], 0);
2170
1.84k
  for (i = 0; i < tab->n_param; 
++i1.32k
) {
2171
1.32k
    if (tab->var[i].is_row)
2172
63
      continue;
2173
1.26k
    col = tab->var[i].index;
2174
1.26k
    isl_int_neg(r_row[off + col], tab->mat->row[row][off + col]);
2175
1.26k
    isl_int_fdiv_r(r_row[off + col], r_row[off + col],
2176
1.26k
        tab->mat->row[row][0]);
2177
1.26k
    isl_int_neg(r_row[off + col], r_row[off + col]);
2178
1.32k
  }
2179
1.69k
  for (i = 0; i < tab->n_div; 
++i1.17k
) {
2180
1.17k
    if (tab->var[tab->n_var - tab->n_div + i].is_row)
2181
31
      continue;
2182
1.14k
    col = tab->var[tab->n_var - tab->n_div + i].index;
2183
1.14k
    isl_int_neg(r_row[off + col], tab->mat->row[row][off + col]);
2184
1.14k
    isl_int_fdiv_r(r_row[off + col], r_row[off + col],
2185
1.14k
        tab->mat->row[row][0]);
2186
1.14k
    isl_int_neg(r_row[off + col], r_row[off + col]);
2187
1.17k
  }
2188
4.01k
  for (i = 0; i < tab->n_col; 
++i3.49k
) {
2189
3.49k
    if (tab->col_var[i] >= 0 &&
2190
3.49k
        
(2.40k
tab->col_var[i] < tab->n_param2.40k
||
2191
2.40k
         
tab->col_var[i] >= tab->n_var - tab->n_div1.14k
))
2192
2.40k
      continue;
2193
1.09k
    isl_int_fdiv_r(r_row[off + i],
2194
1.09k
      tab->mat->row[row][off + i], tab->mat->row[row][0]);
2195
1.09k
  }
2196
517
  if (tab->var[tab->n_var - tab->n_div + d].is_row) {
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
516
  } else {
2209
516
    col = tab->var[tab->n_var - tab->n_div + d].index;
2210
516
    isl_int_set(r_row[off + col], tab->mat->row[row][0]);
2211
516
  }
2212
517
2213
517
  tab->con[r].is_nonneg = 1;
2214
517
  if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r]) < 0)
2215
0
    return -1;
2216
517
  if (tab->row_sign)
2217
517
    tab->row_sign[tab->con[r].index] = isl_tab_row_neg;
2218
517
2219
517
  row = tab->con[r].index;
2220
517
2221
517
  if (d >= n && 
context->op->detect_equalities(context, tab) < 0515
)
2222
0
    return -1;
2223
517
2224
517
  return row;
2225
517
}
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
7.36k
{
2242
7.36k
  int i;
2243
7.36k
  struct isl_tab *tab;
2244
7.36k
  unsigned n_var;
2245
7.36k
  unsigned o_var;
2246
7.36k
2247
7.36k
  tab = isl_tab_alloc(bmap->ctx, 2 * bmap->n_eq + bmap->n_ineq + 1,
2248
7.36k
          isl_basic_map_total_dim(bmap), M);
2249
7.36k
  if (!tab)
2250
0
    return NULL;
2251
7.36k
2252
7.36k
  tab->rational = ISL_F_ISSET(bmap, ISL_BASIC_MAP_RATIONAL);
2253
7.36k
  if (dom) {
2254
6.92k
    tab->n_param = isl_basic_set_total_dim(dom) - dom->n_div;
2255
6.92k
    tab->n_div = dom->n_div;
2256
6.92k
    tab->row_sign = isl_calloc_array(bmap->ctx,
2257
6.92k
          enum isl_tab_row_sign, tab->mat->n_row);
2258
6.92k
    if (tab->mat->n_row && !tab->row_sign)
2259
0
      goto error;
2260
7.36k
  }
2261
7.36k
  if (ISL_F_ISSET(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
7.36k
2267
29.6k
  
for (i = tab->n_param; 7.36k
i < tab->n_var - tab->n_div;
++i22.2k
) {
2268
22.2k
    tab->var[i].is_nonneg = 1;
2269
22.2k
    tab->var[i].frozen = 1;
2270
22.2k
  }
2271
7.36k
  o_var = 1 + tab->n_param;
2272
7.36k
  n_var = tab->n_var - tab->n_param - tab->n_div;
2273
30.3k
  for (i = 0; i < bmap->n_eq; 
++i23.0k
) {
2274
23.0k
    if (max)
2275
18.0k
      isl_seq_neg(bmap->eq[i] + o_var,
2276
18.0k
            bmap->eq[i] + o_var, n_var);
2277
23.0k
    tab = add_lexmin_valid_eq(tab, bmap->eq[i]);
2278
23.0k
    if (max)
2279
18.0k
      isl_seq_neg(bmap->eq[i] + o_var,
2280
18.0k
            bmap->eq[i] + o_var, n_var);
2281
23.0k
    if (!tab || tab->empty)
2282
0
      return tab;
2283
23.0k
  }
2284
7.36k
  if (bmap->n_eq && 
restore_lexmin(tab) < 06.35k
)
2285
0
    goto error;
2286
29.4k
  
for (i = 0; 7.36k
i < bmap->n_ineq;
++i22.1k
) {
2287
22.1k
    if (max)
2288
13.9k
      isl_seq_neg(bmap->ineq[i] + o_var,
2289
13.9k
            bmap->ineq[i] + o_var, n_var);
2290
22.1k
    tab = add_lexmin_ineq(tab, bmap->ineq[i]);
2291
22.1k
    if (max)
2292
13.9k
      isl_seq_neg(bmap->ineq[i] + o_var,
2293
13.9k
            bmap->ineq[i] + o_var, n_var);
2294
22.1k
    if (!tab || tab->empty)
2295
0
      return tab;
2296
22.1k
  }
2297
7.36k
  return tab;
2298
0
error:
2299
0
  isl_tab_free(tab);
2300
0
  return NULL;
2301
7.36k
}
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
341
{
2328
341
  struct isl_tab_undo *snap;
2329
341
  int split;
2330
341
  int row;
2331
341
  int best = -1;
2332
341
  int best_r = 0;
2333
341
2334
341
  if (isl_tab_extend_cons(context_tab, 2) < 0)
2335
0
    return -1;
2336
341
2337
341
  snap = isl_tab_snap(context_tab);
2338
341
2339
2.72k
  for (split = tab->n_redundant; split < tab->n_row; 
++split2.38k
) {
2340
2.38k
    struct isl_tab_undo *snap2;
2341
2.38k
    struct isl_vec *ineq = NULL;
2342
2.38k
    int r = 0;
2343
2.38k
    int ok;
2344
2.38k
2345
2.38k
    if (!isl_tab_var_from_row(tab, split)->is_nonneg)
2346
0
      continue;
2347
2.38k
    if (tab->row_sign[split] != isl_tab_row_any)
2348
1.56k
      continue;
2349
817
2350
817
    ineq = get_row_parameter_ineq(tab, split);
2351
817
    if (!ineq)
2352
0
      return -1;
2353
817
    ok = isl_tab_add_ineq(context_tab, ineq->el) >= 0;
2354
817
    isl_vec_free(ineq);
2355
817
    if (!ok)
2356
0
      return -1;
2357
817
2358
817
    snap2 = isl_tab_snap(context_tab);
2359
817
2360
6.85k
    for (row = tab->n_redundant; row < tab->n_row; 
++row6.04k
) {
2361
6.04k
      struct isl_tab_var *var;
2362
6.04k
2363
6.04k
      if (row == split)
2364
817
        continue;
2365
5.22k
      if (!isl_tab_var_from_row(tab, row)->is_nonneg)
2366
0
        continue;
2367
5.22k
      if (tab->row_sign[row] != isl_tab_row_any)
2368
3.82k
        continue;
2369
1.40k
2370
1.40k
      ineq = get_row_parameter_ineq(tab, row);
2371
1.40k
      if (!ineq)
2372
0
        return -1;
2373
1.40k
      ok = isl_tab_add_ineq(context_tab, ineq->el) >= 0;
2374
1.40k
      isl_vec_free(ineq);
2375
1.40k
      if (!ok)
2376
0
        return -1;
2377
1.40k
      var = &context_tab->con[context_tab->n_con - 1];
2378
1.40k
      if (!context_tab->empty &&
2379
1.40k
          !isl_tab_min_at_most_neg_one(context_tab, var))
2380
26
        r++;
2381
1.40k
      if (isl_tab_rollback(context_tab, snap2) < 0)
2382
0
        return -1;
2383
6.04k
    }
2384
817
    if (best == -1 || 
r > best_r476
) {
2385
346
      best = split;
2386
346
      best_r = r;
2387
346
    }
2388
817
    if (isl_tab_rollback(context_tab, snap) < 0)
2389
0
      return -1;
2390
2.38k
  }
2391
341
2392
341
  return best;
2393
341
}
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 (add_lexmin_eq(clex->tab, eq) < 0)
2417
0
    goto error;
2418
0
  if (check) {
2419
0
    int v = tab_has_valid_sample(clex->tab, eq, 1);
2420
0
    if (v < 0)
2421
0
      goto error;
2422
0
    if (!v)
2423
0
      clex->tab = check_integer_feasible(clex->tab);
2424
0
  }
2425
0
  if (update)
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 (check) {
2441
0
    int v = tab_has_valid_sample(clex->tab, ineq, 0);
2442
0
    if (v < 0)
2443
0
      goto error;
2444
0
    if (!v)
2445
0
      clex->tab = check_integer_feasible(clex->tab);
2446
0
  }
2447
0
  if (update)
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_ok : isl_stat_error;
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
12.1k
{
2468
12.1k
  int i;
2469
12.1k
  int sgn;
2470
12.1k
  isl_int tmp;
2471
12.1k
  enum isl_tab_row_sign res = isl_tab_row_unknown;
2472
12.1k
2473
12.1k
  isl_assert(tab->mat->ctx, tab->samples, return isl_tab_row_unknown);
2474
12.1k
  isl_assert(tab->mat->ctx, tab->samples->n_col == 1 + tab->n_var,
2475
12.1k
      return isl_tab_row_unknown);
2476
12.1k
2477
12.1k
  isl_int_init(tmp);
2478
28.6k
  for (i = tab->n_outside; i < tab->n_sample; 
++i16.5k
) {
2479
17.1k
    isl_seq_inner_product(tab->samples->row[i], ineq,
2480
17.1k
          1 + tab->n_var, &tmp);
2481
17.1k
    sgn = isl_int_sgn(tmp);
2482
17.1k
    if (sgn > 0 || 
(7.69k
sgn == 07.69k
&&
strict4.69k
)) {
2483
13.6k
      if (res == isl_tab_row_unknown)
2484
9.64k
        res = isl_tab_row_pos;
2485
13.6k
      if (res == isl_tab_row_neg)
2486
253
        res = isl_tab_row_any;
2487
13.6k
    }
2488
17.1k
    if (sgn < 0) {
2489
3.00k
      if (res == isl_tab_row_unknown)
2490
2.40k
        res = isl_tab_row_neg;
2491
3.00k
      if (res == isl_tab_row_pos)
2492
300
        res = isl_tab_row_any;
2493
3.00k
    }
2494
17.1k
    if (res == isl_tab_row_any)
2495
553
      break;
2496
17.1k
  }
2497
12.1k
  isl_int_clear(tmp);
2498
12.1k
2499
12.1k
  return res;
2500
12.1k
}
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 (isl_tab_extend_cons(clex->tab, 1) < 0)
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 (clex->tab->M)
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 && isl_tab_rollback(clex->tab, snap) < 0)
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 (isl_tab_save_samples(clex->tab) < 0)
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) < 0) {
2616
0
    isl_tab_free(clex->tab);
2617
0
    clex->tab = NULL;
2618
0
  }
2619
0
}
2620
2621
static void context_lex_discard(void *save)
2622
0
{
2623
0
}
2624
2625
static int context_lex_is_ok(struct isl_context *context)
2626
0
{
2627
0
  struct isl_context_lex *clex = (struct isl_context_lex *)context;
2628
0
  return !!clex->tab;
2629
0
}
2630
2631
/* For each variable in the context tableau, check if the variable can
2632
 * only attain non-negative values.  If so, mark the parameter as non-negative
2633
 * in the main tableau.  This allows for a more direct identification of some
2634
 * cases of violated constraints.
2635
 */
2636
static struct isl_tab *tab_detect_nonnegative_parameters(struct isl_tab *tab,
2637
  struct isl_tab *context_tab)
2638
6.92k
{
2639
6.92k
  int i;
2640
6.92k
  struct isl_tab_undo *snap;
2641
6.92k
  struct isl_vec *ineq = NULL;
2642
6.92k
  struct isl_tab_var *var;
2643
6.92k
  int n;
2644
6.92k
2645
6.92k
  if (context_tab->n_var == 0)
2646
1.45k
    return tab;
2647
5.46k
2648
5.46k
  ineq = isl_vec_alloc(tab->mat->ctx, 1 + context_tab->n_var);
2649
5.46k
  if (!ineq)
2650
0
    goto error;
2651
5.46k
2652
5.46k
  if (isl_tab_extend_cons(context_tab, 1) < 0)
2653
0
    goto error;
2654
5.46k
2655
5.46k
  snap = isl_tab_snap(context_tab);
2656
5.46k
2657
5.46k
  n = 0;
2658
5.46k
  isl_seq_clr(ineq->el, ineq->size);
2659
26.7k
  for (i = 0; i < context_tab->n_var; 
++i21.2k
) {
2660
21.2k
    isl_int_set_si(ineq->el[1 + i], 1);
2661
21.2k
    if (isl_tab_add_ineq(context_tab, ineq->el) < 0)
2662
0
      goto error;
2663
21.2k
    var = &context_tab->con[context_tab->n_con - 1];
2664
21.2k
    if (!context_tab->empty &&
2665
21.2k
        
!isl_tab_min_at_most_neg_one(context_tab, var)20.9k
) {
2666
15.7k
      int j = i;
2667
15.7k
      if (i >= tab->n_param)
2668
100
        j = i - tab->n_param + tab->n_var - tab->n_div;
2669
15.7k
      tab->var[j].is_nonneg = 1;
2670
15.7k
      n++;
2671
15.7k
    }
2672
21.2k
    isl_int_set_si(ineq->el[1 + i], 0);
2673
21.2k
    if (isl_tab_rollback(context_tab, snap) < 0)
2674
0
      goto error;
2675
21.2k
  }
2676
5.46k
2677
5.46k
  if (context_tab->M && 
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.46k
2682
5.46k
  isl_vec_free(ineq);
2683
5.46k
  return tab;
2684
0
error:
2685
0
  isl_vec_free(ineq);
2686
0
  isl_tab_free(tab);
2687
0
  return NULL;
2688
6.92k
}
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 = isl_alloc_type(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.92k
{
2812
6.92k
  struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
2813
6.92k
  if (!tab)
2814
0
    return NULL;
2815
6.92k
  return tab_detect_nonnegative_parameters(tab, cgbr->tab);
2816
6.92k
}
2817
2818
static struct isl_basic_set *context_gbr_peek_basic_set(
2819
  struct isl_context *context)
2820
24.9k
{
2821
24.9k
  struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
2822
24.9k
  if (!cgbr->tab)
2823
0
    return NULL;
2824
24.9k
  return isl_tab_peek_bset(cgbr->tab);
2825
24.9k
}
2826
2827
static struct isl_tab *context_gbr_peek_tab(struct isl_context *context)
2828
32.5k
{
2829
32.5k
  struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
2830
32.5k
  return cgbr->tab;
2831
32.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
169
{
2841
169
  int i, j;
2842
169
  struct isl_vec *cst;
2843
169
  struct isl_basic_set *bset = isl_tab_peek_bset(cgbr->tab);
2844
169
  unsigned dim = isl_basic_set_total_dim(bset);
2845
169
2846
169
  cst = isl_vec_alloc(cgbr->tab->mat->ctx, bset->n_ineq);
2847
169
  if (!cst)
2848
0
    return;
2849
169
2850
1.59k
  
for (i = 0; 169
i < bset->n_ineq;
++i1.42k
) {
2851
1.42k
    isl_int_set(cst->el[i], bset->ineq[i][0]);
2852
10.6k
    for (j = 0; j < dim; 
++j9.21k
) {
2853
9.21k
      if (!isl_int_is_neg(bset->ineq[i][1 + j]))
2854
9.21k
        
continue7.90k
;
2855
1.30k
      isl_int_add(bset->ineq[i][0], bset->ineq[i][0],
2856
1.30k
            bset->ineq[i][1 + j]);
2857
1.30k
    }
2858
1.42k
  }
2859
169
2860
169
  cgbr->shifted = isl_tab_from_basic_set(bset, 0);
2861
169
2862
1.59k
  for (i = 0; i < bset->n_ineq; 
++i1.42k
)
2863
1.42k
    isl_int_set(bset->ineq[i][0], cst->el[i]);
2864
169
2865
169
  isl_vec_free(cst);
2866
169
}
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
210
{
2874
210
  struct isl_vec *sample;
2875
210
2876
210
  if (!cgbr->shifted)
2877
169
    gbr_init_shifted(cgbr);
2878
210
  if (!cgbr->shifted)
2879
0
    return NULL;
2880
210
  if (cgbr->shifted->empty)
2881
161
    return isl_vec_alloc(cgbr->tab->mat->ctx, 0);
2882
49
2883
49
  sample = isl_tab_get_sample_value(cgbr->shifted);
2884
49
  sample = isl_vec_ceil(sample);
2885
49
2886
49
  return sample;
2887
49
}
2888
2889
static __isl_give isl_basic_set *drop_constant_terms(
2890
  __isl_take isl_basic_set *bset)
2891
811
{
2892
811
  int i;
2893
811
2894
811
  if (!bset)
2895
0
    return NULL;
2896
811
2897
2.07k
  
for (i = 0; 811
i < bset->n_eq;
++i1.26k
)
2898
1.26k
    isl_int_set_si(bset->eq[i][0], 0);
2899
811
2900
6.28k
  for (i = 0; i < bset->n_ineq; 
++i5.47k
)
2901
5.47k
    isl_int_set_si(bset->ineq[i][0], 0);
2902
811
2903
811
  return bset;
2904
811
}
2905
2906
static int use_shifted(struct isl_context_gbr *cgbr)
2907
1.04k
{
2908
1.04k
  if (!cgbr->tab)
2909
0
    return 0;
2910
1.04k
  return cgbr->tab->bmap->n_eq == 0 && 
cgbr->tab->bmap->n_div == 0408
;
2911
1.04k
}
2912
2913
static struct isl_vec *gbr_get_sample(struct isl_context_gbr *cgbr)
2914
11.4k
{
2915
11.4k
  struct isl_basic_set *bset;
2916
11.4k
  struct isl_basic_set *cone;
2917
11.4k
2918
11.4k
  if (isl_tab_sample_is_integer(cgbr->tab))
2919
10.4k
    return isl_tab_get_sample_value(cgbr->tab);
2920
1.03k
2921
1.03k
  if (use_shifted(cgbr)) {
2922
210
    struct isl_vec *sample;
2923
210
2924
210
    sample = gbr_get_shifted_sample(cgbr);
2925
210
    if (!sample || sample->size > 0)
2926
49
      return sample;
2927
161
2928
161
    isl_vec_free(sample);
2929
161
  }
2930
1.03k
2931
1.03k
  
if (990
!cgbr->cone990
) {
2932
751
    bset = isl_tab_peek_bset(cgbr->tab);
2933
751
    cgbr->cone = isl_tab_from_recession_cone(bset, 0);
2934
751
    if (!cgbr->cone)
2935
0
      return NULL;
2936
751
    if (isl_tab_track_bset(cgbr->cone,
2937
751
          isl_basic_set_copy(bset)) < 0)
2938
0
      return NULL;
2939
990
  }
2940
990
  if (isl_tab_detect_implicit_equalities(cgbr->cone) < 0)
2941
0
    return NULL;
2942
990
2943
990
  if (cgbr->cone->n_dead == cgbr->cone->n_col) {
2944
179
    struct isl_vec *sample;
2945
179
    struct isl_tab_undo *snap;
2946
179
2947
179
    if (cgbr->tab->basis) {
2948
49
      if (cgbr->tab->basis->n_col != 1 + cgbr->tab->n_var) {
2949
21
        isl_mat_free(cgbr->tab->basis);
2950
21
        cgbr->tab->basis = NULL;
2951
21
      }
2952
49
      cgbr->tab->n_zero = 0;
2953
49
      cgbr->tab->n_unbounded = 0;
2954
49
    }
2955
179
2956
179
    snap = isl_tab_snap(cgbr->tab);
2957
179
2958
179
    sample = isl_tab_sample(cgbr->tab);
2959
179
2960
179
    if (!sample || isl_tab_rollback(cgbr->tab, snap) < 0) {
2961
0
      isl_vec_free(sample);
2962
0
      return NULL;
2963
0
    }
2964
179
2965
179
    return sample;
2966
179
  }
2967
811
2968
811
  cone = isl_basic_set_dup(isl_tab_peek_bset(cgbr->cone));
2969
811
  cone = drop_constant_terms(cone);
2970
811
  cone = isl_basic_set_update_from_tab(cone, cgbr->cone);
2971
811
  cone = isl_basic_set_underlying_set(cone);
2972
811
  cone = isl_basic_set_gauss(cone, NULL);
2973
811
2974
811
  bset = isl_basic_set_dup(isl_tab_peek_bset(cgbr->tab));
2975
811
  bset = isl_basic_set_update_from_tab(bset, cgbr->tab);
2976
811
  bset = isl_basic_set_underlying_set(bset);
2977
811
  bset = isl_basic_set_gauss(bset, NULL);
2978
811
2979
811
  return isl_basic_set_sample_with_cone(bset, cone);
2980
811
}
2981
2982
static void check_gbr_integer_feasible(struct isl_context_gbr *cgbr)
2983
38.2k
{
2984
38.2k
  struct isl_vec *sample;
2985
38.2k
2986
38.2k
  if (!cgbr->tab)
2987
0
    return;
2988
38.2k
2989
38.2k
  if (cgbr->tab->empty)
2990
26.7k
    return;
2991
11.4k
2992
11.4k
  sample = gbr_get_sample(cgbr);
2993
11.4k
  if (!sample)
2994
0
    goto error;
2995
11.4k
2996
11.4k
  if (sample->size == 0) {
2997
252
    isl_vec_free(sample);
2998
252
    if (isl_tab_mark_empty(cgbr->tab) < 0)
2999
0
      goto error;
3000
252
    return;
3001
252
  }
3002
11.2k
3003
11.2k
  if (isl_tab_add_sample(cgbr->tab, sample) < 0)
3004
0
    goto error;
3005
11.2k
3006
11.2k
  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.30k
{
3014
9.30k
  if (!tab)
3015
0
    return NULL;
3016
9.30k
3017
9.30k
  if (isl_tab_extend_cons(tab, 2) < 0)
3018
0
    goto error;
3019
9.30k
3020
9.30k
  if (isl_tab_add_eq(tab, eq) < 0)
3021
0
    goto error;
3022
9.30k
3023
9.30k
  return tab;
3024
0
error:
3025
0
  isl_tab_free(tab);
3026
0
  return NULL;
3027
9.30k
}
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.30k
{
3041
9.30k
  struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3042
9.30k
3043
9.30k
  cgbr->tab = add_gbr_eq(cgbr->tab, eq);
3044
9.30k
3045
9.30k
  if (cgbr->shifted && 
!cgbr->shifted->empty58
&&
use_shifted(cgbr)1
) {
3046
1
    if (isl_tab_mark_empty(cgbr->shifted) < 0)
3047
0
      goto error;
3048
9.30k
  }
3049
9.30k
3050
9.30k
  if (cgbr->cone && 
cgbr->cone->n_col != cgbr->cone->n_dead469
) {
3051
447
    if (isl_tab_extend_cons(cgbr->cone, 2) < 0)
3052
0
      goto error;
3053
447
    if (isl_tab_add_eq(cgbr->cone, eq) < 0)
3054
0
      goto error;
3055
9.30k
  }
3056
9.30k
3057
9.30k
  if (check) {
3058
9.30k
    int v = tab_has_valid_sample(cgbr->tab, eq, 1);
3059
9.30k
    if (v < 0)
3060
0
      goto error;
3061
9.30k
    if (!v)
3062
106
      check_gbr_integer_feasible(cgbr);
3063
9.30k
  }
3064
9.30k
  if (update)
3065
9.30k
    cgbr->tab = check_samples(cgbr->tab, eq, 1);
3066
9.30k
  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
37.4k
{
3074
37.4k
  if (!cgbr->tab)
3075
0
    return;
3076
37.4k
3077
37.4k
  if (isl_tab_extend_cons(cgbr->tab, 1) < 0)
3078
0
    goto error;
3079
37.4k
3080
37.4k
  if (isl_tab_add_ineq(cgbr->tab, ineq) < 0)
3081
0
    goto error;
3082
37.4k
3083
37.4k
  if (cgbr->shifted && 
!cgbr->shifted->empty503
&&
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 < dim;
++i8
) {
3092
8
      if (!isl_int_is_neg(ineq[1 + i]))
3093
8
        
continue6
;
3094
2
      isl_int_add(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 < dim;
++i8
) {
3101
8
      if (!isl_int_is_neg(ineq[1 + i]))
3102
8
        
continue6
;
3103
2
      isl_int_sub(ineq[0], ineq[0], ineq[1 + i]);
3104
2
    }
3105
2
  }
3106
37.4k
3107
37.4k
  if (cgbr->cone && 
cgbr->cone->n_col != cgbr->cone->n_dead4.40k
) {
3108
3.18k
    if (isl_tab_extend_cons(cgbr->cone, 1) < 0)
3109
0
      goto error;
3110
3.18k
    if (isl_tab_add_ineq(cgbr->cone, ineq) < 0)
3111
0
      goto error;
3112
37.4k
  }
3113
37.4k
3114
37.4k
  return;
3115
0
error:
3116
0
  isl_tab_free(cgbr->tab);
3117
0
  cgbr->tab = NULL;
3118
0
}
3119
3120
static void context_gbr_add_ineq(struct isl_context *context, isl_int *ineq,
3121
    int check, int update)
3122
25.8k
{
3123
25.8k
  struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3124
25.8k
3125
25.8k
  add_gbr_ineq(cgbr, ineq);
3126
25.8k
  if (!cgbr->tab)
3127
0
    return;
3128
25.8k
3129
25.8k
  if (check) {
3130
18.7k
    int v = tab_has_valid_sample(cgbr->tab, ineq, 0);
3131
18.7k
    if (v < 0)
3132
0
      goto error;
3133
18.7k
    if (!v)
3134
18.5k
      check_gbr_integer_feasible(cgbr);
3135
18.7k
  }
3136
25.8k
  if (update)
3137
6.00k
    cgbr->tab = check_samples(cgbr->tab, ineq, 0);
3138
25.8k
  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.16k
{
3146
1.16k
  struct isl_context *context = (struct isl_context *)user;
3147
1.16k
  context_gbr_add_ineq(context, ineq, 0, 0);
3148
1.16k
  return context->op->is_ok(context) ? isl_stat_ok : 
isl_stat_error0
;
3149
1.16k
}
3150
3151
static enum isl_tab_row_sign context_gbr_ineq_sign(struct isl_context *context,
3152
      isl_int *ineq, int strict)
3153
12.1k
{
3154
12.1k
  struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3155
12.1k
  return tab_ineq_sign(cgbr->tab, ineq, strict);
3156
12.1k
}
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
11.5k
{
3163
11.5k
  struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3164
11.5k
  struct isl_tab_undo *snap;
3165
11.5k
  struct isl_tab_undo *shifted_snap = NULL;
3166
11.5k
  struct isl_tab_undo *cone_snap = NULL;
3167
11.5k
  int feasible;
3168
11.5k
3169
11.5k
  if (!cgbr->tab)
3170
0
    return -1;
3171
11.5k
3172
11.5k
  if (isl_tab_extend_cons(cgbr->tab, 1) < 0)
3173
0
    return -1;
3174
11.5k
3175
11.5k
  snap = isl_tab_snap(cgbr->tab);
3176
11.5k
  if (cgbr->shifted)
3177
343
    shifted_snap = isl_tab_snap(cgbr->shifted);
3178
11.5k
  if (cgbr->cone)
3179
2.21k
    cone_snap = isl_tab_snap(cgbr->cone);
3180
11.5k
  add_gbr_ineq(cgbr, ineq);
3181
11.5k
  check_gbr_integer_feasible(cgbr);
3182
11.5k
  if (!cgbr->tab)
3183
0
    return -1;
3184
11.5k
  feasible = !cgbr->tab->empty;
3185
11.5k
  if (isl_tab_rollback(cgbr->tab, snap) < 0)
3186
0
    return -1;
3187
11.5k
  if (shifted_snap) {
3188
343
    if (isl_tab_rollback(cgbr->shifted, shifted_snap))
3189
0
      return -1;
3190
11.2k
  } else if (cgbr->shifted) {
3191
114
    isl_tab_free(cgbr->shifted);
3192
114
    cgbr->shifted = NULL;
3193
114
  }
3194
11.5k
  if (cone_snap) {
3195
2.21k
    if (isl_tab_rollback(cgbr->cone, cone_snap))
3196
0
      return -1;
3197
9.36k
  } else if (cgbr->cone) {
3198
283
    isl_tab_free(cgbr->cone);
3199
283
    cgbr->cone = NULL;
3200
283
  }
3201
11.5k
3202
11.5k
  return feasible;
3203
11.5k
}
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
166
{
3212
166
  int i;
3213
166
  int col;
3214
166
3215
166
  if (tab->n_var == 0)
3216
0
    return -1;
3217
166
3218
457
  
for (i = tab->n_var - 1; 166
i >= 0;
--i291
) {
3219
457
    if (i >= tab->n_param && 
i < tab->n_var - tab->n_div440
)
3220
31
      continue;
3221
426
    if (tab->var[i].is_row)
3222
73
      continue;
3223
353
    col = tab->var[i].index;
3224
353
    if (!isl_int_is_zero(p[col]))
3225
353
      
return col166
;
3226
457
  }
3227
166
3228
166
  
return -10
;
3229
166
}
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
104
{
3251
104
  int i;
3252
104
  struct isl_vec *eq = NULL;
3253
104
  isl_bool needs_undo;
3254
104
3255
104
  needs_undo = isl_tab_need_undo(tab);
3256
104
  if (needs_undo < 0)
3257
0
    goto error;
3258
104
  eq = isl_vec_alloc(tab->mat->ctx, 1 + tab->n_var);
3259
104
  if (!eq)
3260
0
    goto error;
3261
104
3262
104
  if (isl_tab_extend_cons(tab, (cgbr->tab->bmap->n_ineq - first)/2) < 0)
3263
0
    goto error;
3264
104
3265
104
  isl_seq_clr(eq->el + 1 + tab->n_param,
3266
104
        tab->n_var - tab->n_param - tab->n_div);
3267
270
  for (i = first; i < cgbr->tab->bmap->n_ineq; 
i += 2166
) {
3268
166
    int j;
3269
166
    int r;
3270
166
    struct isl_tab_undo *snap;
3271
166
    snap = isl_tab_snap(tab);
3272
166
3273
166
    isl_seq_cpy(eq->el, cgbr->tab->bmap->ineq[i], 1 + tab->n_param);
3274
166
    isl_seq_cpy(eq->el + 1 + tab->n_var - tab->n_div,
3275
166
          cgbr->tab->bmap->ineq[i] + 1 + tab->n_param,
3276
166
          tab->n_div);
3277
166
3278
166
    r = isl_tab_add_row(tab, eq->el);
3279
166
    if (r < 0)
3280
0
      goto error;
3281
166
    r = tab->con[r].index;
3282
166
    j = last_non_zero_var_col(tab, tab->mat->row[r] + 2 + tab->M);
3283
166
    if (j < 0 || j < tab->n_dead ||
3284
166
        !isl_int_is_one(tab->mat->row[r][0]) ||
3285
166
        (!isl_int_is_one(tab->mat->row[r][2 + tab->M + j]) &&
3286
166
         
!102
isl_int_is_negone102
(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
113
    if (isl_tab_pivot(tab, r, j) < 0)
3292
0
      goto error;
3293
113
    if (isl_tab_kill_col(tab, j) < 0)
3294
0
      goto error;
3295
113
3296
113
    if (restore_lexmin(tab) < 0)
3297
0
      goto error;
3298
166
  }
3299
104
3300
104
  if (!needs_undo)
3301
104
    isl_tab_clear_undo(tab);
3302
104
  isl_vec_free(eq);
3303
104
3304
104
  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
104
}
3311
3312
static int context_gbr_detect_equalities(struct isl_context *context,
3313
  struct isl_tab *tab)
3314
515
{
3315
515
  struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3316
515
  unsigned n_ineq;
3317
515
3318
515
  if (!cgbr->cone) {
3319
232
    struct isl_basic_set *bset = isl_tab_peek_bset(cgbr->tab);
3320
232
    cgbr->cone = isl_tab_from_recession_cone(bset, 0);
3321
232
    if (!cgbr->cone)
3322
0
      goto error;
3323
232
    if (isl_tab_track_bset(cgbr->cone,
3324
232
          isl_basic_set_copy(bset)) < 0)
3325
0
      goto error;
3326
515
  }
3327
515
  if (isl_tab_detect_implicit_equalities(cgbr->cone) < 0)
3328
0
    goto error;
3329
515
3330
515
  n_ineq = cgbr->tab->bmap->n_ineq;
3331
515
  cgbr->tab = isl_tab_detect_equalities(cgbr->tab, cgbr->cone);
3332
515
  if (!cgbr->tab)
3333
0
    return -1;
3334
515
  if (cgbr->tab->bmap->n_ineq > n_ineq &&
3335
515
      
propagate_equalities(cgbr, tab, n_ineq) < 0104
)
3336
0
    return -1;
3337
515
3338
515
  return 0;
3339
0
error:
3340
0
  isl_tab_free(cgbr->tab);
3341
0
  cgbr->tab = NULL;
3342
0
  return -1;
3343
515
}
3344
3345
static int context_gbr_get_div(struct isl_context *context, struct isl_tab *tab,
3346
    struct isl_vec *div)
3347
607
{
3348
607
  return get_div(tab, context, div);
3349
607
}
3350
3351
static isl_bool context_gbr_insert_div(struct isl_context *context, int pos,
3352
  __isl_keep isl_vec *div)
3353
580
{
3354
580
  struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3355
580
  if (cgbr->cone) {
3356
290
    int r, n_div, o_div;
3357
290
3358
290
    n_div = isl_basic_map_dim(cgbr->cone->bmap, isl_dim_div);
3359
290
    o_div = cgbr->cone->n_var - n_div;
3360
290
3361
290
    if (isl_tab_extend_cons(cgbr->cone, 3) < 0)
3362
0
      return isl_bool_error;
3363
290
    if (isl_tab_extend_vars(cgbr->cone, 1) < 0)
3364
0
      return isl_bool_error;
3365
290
    if ((r = isl_tab_insert_var(cgbr->cone, pos)) <0)
3366
0
      return isl_bool_error;
3367
290
3368
290
    cgbr->cone->bmap = isl_basic_map_insert_div(cgbr->cone->bmap,
3369
290
                r - o_div, div);
3370
290
    if (!cgbr->cone->bmap)
3371
0
      return isl_bool_error;
3372
290
    if (isl_tab_push_var(cgbr->cone, isl_tab_undo_bmap_div,
3373
290
            &cgbr->cone->var[r]) < 0)
3374
0
      return isl_bool_error;
3375
580
  }
3376
580
  return context_tab_insert_div(cgbr->tab, pos, div,
3377
580
          context_gbr_add_ineq_wrap, context);
3378
580
}
3379
3380
static int context_gbr_best_split(struct isl_context *context,
3381
    struct isl_tab *tab)
3382
341
{
3383
341
  struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3384
341
  struct isl_tab_undo *snap;
3385
341
  int r;
3386
341
3387
341
  snap = isl_tab_snap(cgbr->tab);
3388
341
  r = best_split(tab, cgbr->tab);
3389
341
3390
341
  if (r >= 0 && isl_tab_rollback(cgbr->tab, snap) < 0)
3391
0
    return -1;
3392
341
3393
341
  return r;
3394
341
}
3395
3396
static int context_gbr_is_empty(struct isl_context *context)
3397
45.0k
{
3398
45.0k
  struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3399
45.0k
  if (!cgbr->tab)
3400
0
    return -1;
3401
45.0k
  return cgbr->tab->empty;
3402
45.0k
}
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
28.5k
{
3412
28.5k
  struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3413
28.5k
  struct isl_gbr_tab_undo *snap;
3414
28.5k
3415
28.5k
  if (!cgbr->tab)
3416
0
    return NULL;
3417
28.5k
3418
28.5k
  snap = isl_alloc_type(cgbr->tab->mat->ctx, struct isl_gbr_tab_undo);
3419
28.5k
  if (!snap)
3420
0
    return NULL;
3421
28.5k
3422
28.5k
  snap->tab_snap = isl_tab_snap(cgbr->tab);
3423
28.5k
  if (isl_tab_save_samples(cgbr->tab) < 0)
3424
0
    goto error;
3425
28.5k
3426
28.5k
  if (cgbr->shifted)
3427
178
    snap->shifted_snap = isl_tab_snap(cgbr->shifted);
3428
28.3k
  else
3429
28.3k
    snap->shifted_snap = NULL;
3430
28.5k
3431
28.5k
  if (cgbr->cone)
3432
1.48k
    snap->cone_snap = isl_tab_snap(cgbr->cone);
3433
27.0k
  else
3434
27.0k
    snap->cone_snap = NULL;
3435
28.5k
3436
28.5k
  return snap;
3437
0
error:
3438
0
  free(snap);
3439
0
  return NULL;
3440
28.5k
}
3441
3442
static void context_gbr_restore(struct isl_context *context, void *save)
3443
23.5k
{
3444
23.5k
  struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3445
23.5k
  struct isl_gbr_tab_undo *snap = (struct isl_gbr_tab_undo *)save;
3446
23.5k
  if (!snap)
3447
0
    goto error;
3448
23.5k
  if (isl_tab_rollback(cgbr->tab, snap->tab_snap) < 0)
3449
0
    goto error;
3450
23.5k
3451
23.5k
  if (snap->shifted_snap) {
3452
125
    if (isl_tab_rollback(cgbr->shifted, snap->shifted_snap) < 0)
3453
0
      goto error;
3454
23.4k
  } else if (cgbr->shifted) {
3455
0
    isl_tab_free(cgbr->shifted);
3456
0
    cgbr->shifted = NULL;
3457
0
  }
3458
23.5k
3459
23.5k
  if (snap->cone_snap) {
3460
1.40k
    if (isl_tab_rollback(cgbr->cone, snap->cone_snap) < 0)
3461
0
      goto error;
3462
22.1k
  } else if (cgbr->cone) {
3463
258
    isl_tab_free(cgbr->cone);
3464
258
    cgbr->cone = NULL;
3465
258
  }
3466
23.5k
3467
23.5k
  free(snap);
3468
23.5k
3469
23.5k
  return;
3470
0
error:
3471
0
  free(snap);
3472
0
  isl_tab_free(cgbr->tab);
3473
0
  cgbr->tab = NULL;
3474
0
}
3475
3476
static void context_gbr_discard(void *save)
3477
5.02k
{
3478
5.02k
  struct isl_gbr_tab_undo *snap = (struct isl_gbr_tab_undo *)save;
3479
5.02k
  free(snap);
3480
5.02k
}
3481
3482
static int context_gbr_is_ok(struct isl_context *context)
3483
1.83k
{
3484
1.83k
  struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3485
1.83k
  return !!cgbr->tab;
3486
1.83k
}
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.98k
{
3498
7.98k
  struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3499
7.98k
  isl_tab_free(cgbr->tab);
3500
7.98k
  isl_tab_free(cgbr->shifted);
3501
7.98k
  isl_tab_free(cgbr->cone);
3502
7.98k
  free(cgbr);
3503
7.98k
3504
7.98k
  return NULL;
3505
7.98k
}
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.98k
{
3530
7.98k
  struct isl_context_gbr *cgbr;
3531
7.98k
3532
7.98k
  if (!dom)
3533
0
    return NULL;
3534
7.98k
3535
7.98k
  cgbr = isl_calloc_type(dom->ctx, struct isl_context_gbr);
3536
7.98k
  if (!cgbr)
3537
0
    return NULL;
3538
7.98k
3539
7.98k
  cgbr->context.op = &isl_context_gbr_op;
3540
7.98k
3541
7.98k
  cgbr->shifted = NULL;
3542
7.98k
  cgbr->cone = NULL;
3543
7.98k
  cgbr->tab = isl_tab_from_basic_set(dom, 1);
3544
7.98k
  cgbr->tab = isl_tab_init_samples(cgbr->tab);
3545
7.98k
  if (!cgbr->tab)
3546
0
    goto error;
3547
7.98k
  check_gbr_integer_feasible(cgbr);
3548
7.98k
3549
7.98k
  return &cgbr->context;
3550
0
error:
3551
0
  cgbr->context.op->free(&cgbr->context);
3552
0
  return NULL;
3553
7.98k
}
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.98k
{
3563
7.98k
  struct isl_context *context;
3564
7.98k
  int first;
3565
7.98k
3566
7.98k
  if (!dom)
3567
0
    return NULL;
3568
7.98k
3569
7.98k
  if (dom->ctx->opt->context == ISL_CONTEXT_LEXMIN)
3570
7.98k
    
context = isl_context_lex_alloc(dom)0
;
3571
7.98k
  else
3572
7.98k
    context = isl_context_gbr_alloc(dom);
3573
7.98k
3574
7.98k
  if (!context)
3575
0
    return NULL;
3576
7.98k
3577
7.98k
  first = isl_basic_set_first_unknown_div(dom);
3578
7.98k
  if (first < 0)
3579
0
    return context->op->free(context);
3580
7.98k
  context->n_unknown = isl_basic_set_dim(dom, isl_dim_div) - first;
3581
7.98k
3582
7.98k
  return context;
3583
7.98k
}
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.98k
{
3595
7.98k
  sol->rational = ISL_F_ISSET(bmap, ISL_BASIC_MAP_RATIONAL);
3596
7.98k
  sol->dec_level.callback.run = &sol_dec_level_wrap;
3597
7.98k
  sol->dec_level.sol = sol;
3598
7.98k
  sol->max = max;
3599
7.98k
  sol->n_out = isl_basic_map_dim(bmap, isl_dim_out);
3600
7.98k
  sol->space = isl_basic_map_get_space(bmap);
3601
7.98k
3602
7.98k
  sol->context = isl_context_alloc(dom);
3603
7.98k
  if (!sol->space || !sol->context)
3604
0
    return isl_stat_error;
3605
7.98k
3606
7.98k
  return isl_stat_ok;
3607
7.98k
}
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.48k
{
3619
3.48k
  struct isl_sol_map *sol_map = NULL;
3620
3.48k
  isl_space *space;
3621
3.48k
3622
3.48k
  if (!bmap)
3623
0
    goto error;
3624
3.48k
3625
3.48k
  sol_map = isl_calloc_type(bmap->ctx, struct isl_sol_map);
3626
3.48k
  if (!sol_map)
3627
0
    goto error;
3628
3.48k
3629
3.48k
  sol_map->sol.free = &sol_map_free;
3630
3.48k
  if (sol_init(&sol_map->sol, bmap, dom, max) < 0)
3631
0
    goto error;
3632
3.48k
  sol_map->sol.add = &sol_map_add_wrap;
3633
3.48k
  sol_map->sol.add_empty = track_empty ? 
&sol_map_add_empty_wrap2.63k
: NULL;
3634
3.48k
  space = isl_space_copy(sol_map->sol.space);
3635
3.48k
  sol_map->map = isl_map_alloc_space(space, 1, ISL_MAP_DISJOINT);
3636
3.48k
  if (!sol_map->map)
3637
0
    goto error;
3638
3.48k
3639
3.48k
  if (track_empty) {
3640
2.63k
    sol_map->empty = isl_set_alloc_space(isl_basic_set_get_space(dom),
3641
2.63k
              1, ISL_SET_DISJOINT);
3642
2.63k
    if (!sol_map->empty)
3643
0
      goto error;
3644
3.48k
  }
3645
3.48k
3646
3.48k
  isl_basic_set_free(dom);
3647
3.48k
  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.48k
}
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
12.1k
{
3659
12.1k
  int j;
3660
12.1k
  unsigned off = 2 + tab->M;
3661
12.1k
3662
52.9k
  for (j = tab->n_dead; j < tab->n_col; 
++j40.7k
) {
3663
42.1k
    if (tab->col_var[j] >= 0 &&
3664
42.1k
        
(29.9k
tab->col_var[j] < tab->n_param29.9k
||
3665
29.9k
        
tab->col_var[j] >= tab->n_var - tab->n_div1.49k
))
3666
29.9k
      continue;
3667
12.2k
3668
12.2k
    if (isl_int_is_pos(tab->mat->row[row][off + j]))
3669
12.2k
      
return 01.38k
;
3670
42.1k
  }
3671
12.1k
3672
12.1k
  
return 110.7k
;
3673
12.1k
}
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
15.0k
{
3683
15.0k
  isl_int gcd;
3684
15.0k
  int strict = 0;
3685
15.0k
3686
15.0k
  isl_int_init(gcd);
3687
15.0k
  isl_seq_gcd(vec->el + 1, vec->size - 1, &gcd);
3688
15.0k
  if (!isl_int_is_one(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
15.0k
  isl_int_clear(gcd);
3694
15.0k
3695
15.0k
  return strict;
3696
15.0k
}
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
64.9k
{
3757
64.9k
  struct isl_vec *ineq = NULL;
3758
64.9k
  enum isl_tab_row_sign res = isl_tab_row_unknown;
3759
64.9k
  int critical;
3760
64.9k
  int strict;
3761
64.9k
  int row2;
3762
64.9k
3763
64.9k
  if (tab->row_sign[row] != isl_tab_row_unknown)
3764
32.5k
    return tab->row_sign[row];
3765
32.4k
  if (is_obviously_nonneg(tab, row))
3766
20.2k
    return isl_tab_row_pos;
3767
111k
  
for (row2 = tab->n_redundant; 12.1k
row2 < tab->n_row;
++row299.2k
) {
3768
99.2k
    if (tab->row_sign[row2] == isl_tab_row_unknown)
3769
32.0k
      continue;
3770
67.1k
    if (identical_parameter_line(tab, row, row2))
3771
43
      return tab->row_sign[row2];
3772
99.2k
  }
3773
12.1k
3774
12.1k
  critical = is_critical(tab, row);
3775
12.1k
3776
12.1k
  ineq = get_row_parameter_ineq(tab, row);
3777
12.1k
  if (!ineq)
3778
0
    goto error;
3779
12.1k
3780
12.1k
  strict = is_strict(ineq);
3781
12.1k
3782
12.1k
  res = sol->context->op->ineq_sign(sol->context, ineq->el,
3783
12.1k
            critical || 
strict1.38k
);
3784
12.1k
3785
12.1k
  if (res == isl_tab_row_unknown || 
res == isl_tab_row_pos12.0k
) {
3786
9.40k
    /* test for negative values */
3787
9.40k
    int feasible;
3788
9.40k
    isl_seq_neg(ineq->el, ineq->el, ineq->size);
3789
9.40k
    isl_int_sub_ui(ineq->el[0], ineq->el[0], 1);
3790
9.40k
3791
9.40k
    feasible = sol->context->op->test_ineq(sol->context, ineq->el);
3792
9.40k
    if (feasible < 0)
3793
0
      goto error;
3794
9.40k
    if (!feasible)
3795
8.67k
      res = isl_tab_row_pos;
3796
738
    else
3797
738
      res = (res == isl_tab_row_unknown) ? 
isl_tab_row_neg18
3798
738
                 : 
isl_tab_row_any720
;
3799
9.40k
    if (res == isl_tab_row_neg) {
3800
18
      isl_seq_neg(ineq->el, ineq->el, ineq->size);
3801
18
      isl_int_sub_ui(ineq->el[0], ineq->el[0], 1);
3802
18
    }
3803
9.40k
  }
3804
12.1k
3805
12.1k
  if (res == isl_tab_row_neg) {
3806
2.17k
    /* test for positive values */
3807
2.17k
    int feasible;
3808
2.17k
    if (!critical && 
!strict88
)
3809
2.17k
      
isl_int_sub_ui88
(ineq->el[0], ineq->el[0], 1);
3810
2.17k
3811
2.17k
    feasible = sol->context->op->test_ineq(sol->context, ineq->el);
3812
2.17k
    if (feasible < 0)
3813
0
      goto error;
3814
2.17k
    if (feasible)
3815
2.15k
      res = isl_tab_row_any;
3816
2.17k
  }
3817
12.1k
3818
12.1k
  isl_vec_free(ineq);
3819
12.1k
  return res;
3820
0
error:
3821
0
  isl_vec_free(ineq);
3822
0
  return isl_tab_row_unknown;
3823
64.9k
}
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.95k
{
3843
2.95k
  void *saved;
3844
2.95k
3845
2.95k
  if (!sol->context)
3846
0
    goto error;
3847
2.95k
  saved = sol->context->op->save(sol->context);
3848
2.95k
3849
2.95k
  tab = isl_tab_dup(tab);
3850
2.95k
  if (!tab)
3851
0
    goto error;
3852
2.95k
3853
2.95k
  sol->context->op->add_ineq(sol->context, ineq, 0, 1);
3854
2.95k
3855
2.95k
  find_solutions(sol, tab);
3856
2.95k
3857
2.95k
  if (!sol->error)
3858
2.95k
    sol->context->op->restore(sol->context, saved);
3859
0
  else
3860
0
    sol->context->op->discard(saved);
3861
2.95k
  return;
3862
0
error:
3863
0
  sol->error = 1;
3864
0
}
3865
3866
/* Record the absence of solutions for those values of the parameters
3867
 * that do not satisfy the given inequality with equality.
3868
 */
3869
static void no_sol_in_strict(struct isl_sol *sol,
3870
  struct isl_tab *tab, struct isl_vec *ineq)
3871
18.6k
{
3872
18.6k
  int empty;
3873
18.6k
  void *saved;
3874
18.6k
3875
18.6k
  if (!sol->context || sol->error)
3876
0
    goto error;
3877
18.6k
  saved = sol->context->op->save(sol->context);
3878
18.6k
3879
18.6k
  isl_int_sub_ui(ineq->el[0], ineq->el[0], 1);
3880
18.6k
3881
18.6k
  sol->context->op->add_ineq(sol->context, ineq->el, 1, 0);
3882
18.6k
  if (!sol->context)
3883
0
    goto error;
3884
18.6k
3885
18.6k
  empty = tab->empty;
3886
18.6k
  tab->empty = 1;
3887
18.6k
  sol_add(sol, tab);
3888
18.6k
  tab->empty = empty;
3889
18.6k
3890
18.6k
  isl_int_add_ui(ineq->el[0], ineq->el[0], 1);
3891
18.6k
3892
18.6k
  sol->context->op->restore(sol->context, saved);
3893
18.6k
  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.95k
{
3903
2.95k
  int row;
3904
2.95k
3905
23.1k
  for (row = tab->n_redundant; row < tab->n_row; 
++row20.1k
) {
3906
20.1k
    if (!isl_tab_var_from_row(tab, row)->is_nonneg)
3907
17
      continue;
3908
20.1k
    if (tab->row_sign[row] == isl_tab_row_any)
3909
3.43k
      tab->row_sign[row] = isl_tab_row_unknown;
3910
20.1k
  }
3911
2.95k
}
3912
3913
/* Compute the lexicographic minimum of the set represented by the main
3914
 * tableau "tab" within the context "sol->context_tab".
3915
 * On entry the sample value of the main tableau is lexicographically
3916
 * less than or equal to this lexicographic minimum.
3917
 * Pivots are performed until a feasible point is found, which is then
3918
 * necessarily equal to the minimum, or until the tableau is found to
3919
 * be infeasible.  Some pivots may need to be performed for only some
3920
 * feasible values of the context tableau.  If so, the context tableau
3921
 * is split into a part where the pivot is needed and a part where it is not.
3922
 *
3923
 * Whenever we enter the main loop, the main tableau is such that no
3924
 * "obvious" pivots need to be performed on it, where "obvious" means
3925
 * that the given row can be seen to be negative without looking at
3926
 * the context tableau.  In particular, for non-parametric problems,
3927
 * no pivots need to be performed on the main tableau.
3928
 * The caller of find_solutions is responsible for making this property
3929
 * hold prior to the first iteration of the loop, while restore_lexmin
3930
 * is called before every other iteration.
3931
 *
3932
 * Inside the main loop, we first examine the signs of the rows of
3933
 * the main tableau within the context of the context tableau.
3934
 * If we find a row that is always non-positive for all values of
3935
 * the parameters satisfying the context tableau and negative for at
3936
 * least one value of the parameters, we perform the appropriate pivot
3937
 * and start over.  An exception is the case where no pivot can be
3938
 * performed on the row.  In this case, we require that the sign of
3939
 * the row is negative for all values of the parameters (rather than just
3940
 * non-positive).  This special case is handled inside row_sign, which
3941
 * will say that the row can have any sign if it determines that it can
3942
 * attain both negative and zero values.
3943
 *
3944
 * If we can't find a row that always requires a pivot, but we can find
3945
 * one or more rows that require a pivot for some values of the parameters
3946
 * (i.e., the row can attain both positive and negative signs), then we split
3947
 * the context tableau into two parts, one where we force the sign to be
3948
 * non-negative and one where we force is to be negative.
3949
 * The non-negative part is handled by a recursive call (through find_in_pos).
3950
 * Upon returning from this call, we continue with the negative part and
3951
 * perform the required pivot.
3952
 *
3953
 * If no such rows can be found, all rows are non-negative and we have
3954
 * found a (rational) feasible point.  If we only wanted a rational point
3955
 * then we are done.
3956
 * Otherwise, we check if all values of the sample point of the tableau
3957
 * are integral for the variables.  If so, we have found the minimal
3958
 * integral point and we are done.
3959
 * If the sample point is not integral, then we need to make a distinction
3960
 * based on whether the constant term is non-integral or the coefficients
3961
 * of the parameters.  Furthermore, in order to decide how to handle
3962
 * the non-integrality, we also need to know whether the coefficients
3963
 * of the other columns in the tableau are integral.  This leads
3964
 * to the following table.  The first two rows do not correspond
3965
 * to a non-integral sample point and are only mentioned for completeness.
3966
 *
3967
 *  constant  parameters  other
3968
 *
3969
 *  int   int   int |
3970
 *  int   int   rat | -> no problem
3971
 *
3972
 *  rat   int   int   -> fail
3973
 *
3974
 *  rat   int   rat   -> cut
3975
 *
3976
 *  int   rat   rat |
3977
 *  rat   rat   rat | -> parametric cut
3978
 *
3979
 *  int   rat   int |
3980
 *  rat   rat   int | -> split context
3981
 *
3982
 * If the parametric constant is completely integral, then there is nothing
3983
 * to be done.  If the constant term is non-integral, but all the other
3984
 * coefficient are integral, then there is nothing that can be done
3985
 * and the tableau has no integral solution.
3986
 * If, on the other hand, one or more of the other columns have rational
3987
 * coefficients, but the parameter coefficients are all integral, then
3988
 * we can perform a regular (non-parametric) cut.
3989
 * Finally, if there is any parameter coefficient that is non-integral,
3990
 * then we need to involve the context tableau.  There are two cases here.
3991
 * If at least one other column has a rational coefficient, then we
3992
 * can perform a parametric cut in the main tableau by adding a new
3993
 * integer division in the context tableau.
3994
 * If all other columns have integral coefficients, then we need to
3995
 * enforce that the rational combination of parameters (c + \sum a_i y_i)/m
3996
 * is always integral.  We do this by introducing an integer division
3997
 * q = floor((c + \sum a_i y_i)/m) and stipulating that its argument should
3998
 * always be integral in the context tableau, i.e., m q = c + \sum a_i y_i.
3999
 * Since q is expressed in the tableau as
4000
 *  c + \sum a_i y_i - m q >= 0
4001
 *  -c - \sum a_i y_i + m q + m - 1 >= 0
4002
 * it is sufficient to add the inequality
4003
 *  -c - \sum a_i y_i + m q >= 0
4004
 * In the part of the context where this inequality does not hold, the
4005
 * main tableau is marked as being empty.
4006
 */
4007
static void find_solutions(struct isl_sol *sol, struct isl_tab *tab)
4008
9.87k
{
4009
9.87k
  struct isl_context *context;
4010
9.87k
  int r;
4011
9.87k
4012
9.87k
  if (!tab || sol->error)
4013
0
    goto error;
4014
9.87k
4015
9.87k
  context = sol->context;
4016
9.87k
4017
9.87k
  if (tab->empty)
4018
0
    goto done;
4019
9.87k
  if (context->op->is_empty(context))
4020
0
    goto done;
4021
9.87k
4022
13.8k
  
for (r = 0; 9.87k
r >= 0 && tab && !tab->empty;
r = restore_lexmin(tab)4.01k
) {
4023
11.1k
    int flags;
4024
11.1k
    int row;
4025
11.1k
    enum isl_tab_row_sign sgn;
4026
11.1k
    int split = -1;
4027
11.1k
    int n_split = 0;
4028
11.1k
4029
76.8k
    for (row = tab->n_redundant; row < tab->n_row; 
++row65.6k
) {
4030
65.7k
      if (!isl_tab_var_from_row(tab, row)->is_nonneg)
4031
707
        continue;
4032
64.9k
      sgn = row_sign(tab, sol, row);
4033
64.9k
      if (!sgn)
4034
0
        goto error;
4035
64.9k
      tab->row_sign[row] = sgn;
4036
64.9k
      if (sgn == isl_tab_row_any)
4037
3.43k
        n_split++;
4038
64.9k
      if (sgn == isl_tab_row_any && 
split == -13.43k
)
4039
2.95k
        split = row;
4040
64.9k
      if (sgn == isl_tab_row_neg)
4041
18
        break;
4042
65.7k
    }
4043
11.1k
    if (row < tab->n_row)
4044
18
      continue;
4045
11.1k
    if (split != -1) {
4046
2.95k
      struct isl_vec *ineq;
4047
2.95k
      if (n_split != 1)
4048
341
        split = context->op->best_split(context, tab);
4049
2.95k
      if (split < 0)
4050
0
        goto error;
4051
2.95k
      ineq = get_row_parameter_ineq(tab, split);
4052
2.95k
      if (!ineq)
4053
0
        goto error;
4054
2.95k
      is_strict(ineq);
4055
2.95k
      reset_any_to_unknown(tab);
4056
2.95k
      tab->row_sign[split] = isl_tab_row_pos;
4057
2.95k
      sol_inc_level(sol);
4058
2.95k
      find_in_pos(sol, tab, ineq->el);
4059
2.95k
      tab->row_sign[split] = isl_tab_row_neg;
4060
2.95k
      isl_seq_neg(ineq->el, ineq->el, ineq->size);
4061
2.95k
      isl_int_sub_ui(ineq->el[0], ineq->el[0], 1);
4062
2.95k
      if (!sol->error)
4063
2.95k
        context->op->add_ineq(context, ineq->el, 0, 1);
4064
2.95k
      isl_vec_free(ineq);
4065
2.95k
      if (sol->error)
4066
0
        goto error;
4067
2.95k
      continue;
4068
2.95k
    }
4069
8.20k
    if (tab->rational)
4070
3
      break;
4071
8.20k
    row = first_non_integer_row(tab, &flags);
4072
8.20k
    if (row < 0)
4073
7.16k
      break;
4074
1.04k
    if (ISL_FL_ISSET(flags, I_PAR)) {
4075
433
      if (ISL_FL_ISSET(flags, I_VAR)) {
4076
0
        if (isl_tab_mark_empty(tab) < 0)
4077
0
          goto error;
4078
0
        break;
4079
0
      }
4080
433
      row = add_cut(tab, row);
4081
607
    } else if (ISL_FL_ISSET(flags, I_VAR)) {
4082
90
      struct isl_vec *div;
4083
90
      struct isl_vec *ineq;
4084
90
      int d;
4085
90
      div = get_row_split_div(tab, row);
4086
90
      if (!div)
4087
0
        goto error;
4088
90
      d = context->op->get_div(context, tab, div);
4089
90
      isl_vec_free(div);
4090
90
      if (d < 0)
4091
0
        goto error;
4092
90
      ineq = ineq_for_div(context->op->peek_basic_set(context), d);
4093
90
      if (!ineq)
4094
0
        goto error;
4095
90
      sol_inc_level(sol);
4096
90
      no_sol_in_strict(sol, tab, ineq);
4097
90
      isl_seq_neg(ineq->el, ineq->el, ineq->size);
4098
90
      context->op->add_ineq(context, ineq->el, 1, 1);
4099
90
      isl_vec_free(ineq);
4100
90
      if (sol->error || !context->op->is_ok(context))
4101
0
        goto error;
4102
90
      tab = set_row_cst_to_div(tab, row, d);
4103
90
      if (context->op->is_empty(context))
4104
0
        break;
4105
517
    } else
4106
517
      row = add_parametric_cut(tab, row, context);
4107
1.04k
    if (row < 0)
4108
0
      goto error;
4109
11.1k
  }
4110
9.87k
  if (r < 0)
4111
0
    goto error;
4112
9.87k
done:
4113
9.87k
  sol_add(sol, tab);
4114
9.87k
  isl_tab_free(tab);
4115
9.87k
  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.92k
{
4130
6.92k
  if (sol->error)
4131
0
    return 0;
4132
6.92k
  if (!sol->partial)
4133
104
    return 0;
4134
6.81k
  if (!sol->partial->next)
4135
4.83k
    return 0;
4136
1.98k
  return sol->partial->level == sol->partial->next->level;
4137
1.98k
}
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.92k
{
4161
6.92k
  int row;
4162
6.92k
  void *saved;
4163
6.92k
4164
6.92k
  if (!tab)
4165
0
    goto error;
4166
6.92k
4167
6.92k
  sol->level = 0;
4168
6.92k
4169
75.9k
  for (row = tab->n_redundant; row < tab->n_row; 
++row69.0k
) {
4170
69.0k
    int p;
4171
69.0k
    struct isl_vec *eq;
4172
69.0k
4173
69.0k
    if (tab->row_var[row] < 0)
4174
14.1k
      continue;
4175
54.8k
    if (tab->row_var[row] >= tab->n_param &&
4176
54.8k
        
tab->row_var[row] < tab->n_var - tab->n_div45.5k
)
4177
45.5k
      continue;
4178
9.30k
    if (tab->row_var[row] < tab->n_param)
4179
9.30k
      p = tab->row_var[row];
4180
0
    else
4181
0
      p = tab->row_var[row]
4182
0
        + tab->n_param - (tab->n_var - tab->n_div);
4183
9.30k
4184
9.30k
    eq = isl_vec_alloc(tab->mat->ctx, 1+tab->n_param+tab->n_div);
4185
9.30k
    if (!eq)
4186
0
      goto error;
4187
9.30k
    get_row_parameter_line(tab, row, eq->el);
4188
9.30k
    isl_int_neg(eq->el[1 + p], tab->mat->row[row][0]);
4189
9.30k
    eq = isl_vec_normalize(eq);
4190
9.30k
4191
9.30k
    sol_inc_level(sol);
4192
9.30k
    no_sol_in_strict(sol, tab, eq);
4193
9.30k
4194
9.30k
    isl_seq_neg(eq->el, eq->el, eq->size);
4195
9.30k
    sol_inc_level(sol);
4196
9.30k
    no_sol_in_strict(sol, tab, eq);
4197
9.30k
    isl_seq_neg(eq->el, eq->el, eq->size);
4198
9.30k
4199
9.30k
    sol->context->op->add_eq(sol->context, eq->el, 1, 1);
4200
9.30k
4201
9.30k
    isl_vec_free(eq);
4202
9.30k
4203
9.30k
    if (isl_tab_mark_redundant(tab, row) < 0)
4204
0
      goto error;
4205
9.30k
4206
9.30k
    if (sol->context->op->is_empty(sol->context))
4207
0
      break;
4208
9.30k
4209
9.30k
    row = tab->n_redundant - 1;
4210
9.30k
  }
4211
6.92k
4212
6.92k
  saved = sol->context->op->save(sol->context);
4213
6.92k
4214
6.92k
  find_solutions(sol, tab);
4215
6.92k
4216
6.92k
  if (sol_has_mergeable_solutions(sol))
4217
1.89k
    sol->context->op->restore(sol->context, saved);
4218
5.02k
  else
4219
5.02k
    sol->context->op->discard(saved);
4220
6.92k
4221
6.92k
  sol->level = 0;
4222
6.92k
  sol_pop(sol);
4223
6.92k
4224
6.92k
  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
632
{
4237
632
  int i;
4238
632
  unsigned b_dim = isl_space_dim(bmap->dim, isl_dim_all);
4239
632
  unsigned d_dim = isl_space_dim(dom->dim, isl_dim_all);
4240
632
4241
632
  if (isl_int_is_zero(dom->div[div][0]))
4242
632
    
return -1142
;
4243
490
  if (isl_seq_first_non_zero(dom->div[div] + 2 + d_dim, dom->n_div) != -1)
4244
2
    return -1;
4245
488
4246
576
  
for (i = 0; 488
i < bmap->n_div;
++i88
) {
4247
144
    if (isl_int_is_zero(bmap->div[i][0]))
4248
144
      
continue12
;
4249
132
    if (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 (isl_seq_eq(bmap->div[i], dom->div[div], 2 + d_dim))
4253
56
      return i;
4254
144
  }
4255
488
  
return -1432
;
4256
632
}
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
274
{
4275
274
  int i;
4276
274
  int common = 0;
4277
274
  int other;
4278
274
4279
590
  for (i = 0; i < dom->n_div; 
++i316
)
4280
316
    if (find_context_div(bmap, dom, i) != -1)
4281
28
      common++;
4282
274
  other = bmap->n_div - common;
4283
274
  if (dom->n_div - common > 0) {
4284
250
    bmap = isl_basic_map_extend_space(bmap, isl_space_copy(bmap->dim),
4285
250
        dom->n_div - common, 0, 0);
4286
250
    if (!bmap)
4287
0
      return NULL;
4288
274
  }
4289
590
  
for (i = 0; 274
i < dom->n_div;
++i316
) {
4290
316
    int pos = find_context_div(bmap, dom, i);
4291
316
    if (pos < 0) {
4292
288
      pos = isl_basic_map_alloc_div(bmap);
4293
288
      if (pos < 0)
4294
0
        goto error;
4295
288
      isl_int_set_si(bmap->div[pos][0], 0);
4296
288
    }
4297
316
    if (pos != other + i)
4298
15
      isl_basic_map_swap_div(bmap, pos, other + i);
4299
316
  }
4300
274
  return bmap;
4301
0
error:
4302
0
  isl_basic_map_free(bmap);
4303
0
  return NULL;
4304
274
}
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
7.23k
{
4324
7.23k
  struct isl_tab *tab;
4325
7.23k
  struct isl_sol *sol = NULL;
4326
7.23k
  struct isl_context *context;
4327
7.23k
4328
7.23k
  if (dom->n_div) {
4329
274
    dom = isl_basic_set_sort_divs(dom);
4330
274
    bmap = align_context_divs(bmap, dom);
4331
274
  }
4332
7.23k
  sol = init(bmap, dom, !!empty, max);
4333
7.23k
  if (!sol)
4334
0
    goto error;
4335
7.23k
4336
7.23k
  context = sol->context;
4337
7.23k
  if (isl_basic_set_plain_is_empty(context->op->peek_basic_set(context)))
4338
12
    /* nothing */;
4339
7.22k
  else if (isl_basic_map_plain_is_empty(bmap)) {
4340
1.05k
    if (sol->add_empty)
4341
1.04k
      sol->add_empty(sol,
4342
1.04k
        isl_basic_set_copy(context->op->peek_basic_set(context)));
4343
6.16k
  } else {
4344
6.16k
    tab = tab_for_lexmin(bmap,
4345
6.16k
            context->op->peek_basic_set(context), 1, max);
4346
6.16k
    tab = context->op->detect_nonnegative_parameters(context, tab);
4347
6.16k
    find_solutions_main(sol, tab);
4348
6.16k
  }
4349
7.23k
  if (sol->error)
4350
0
    goto error;
4351
7.23k
4352
7.23k
  isl_basic_map_free(bmap);
4353
7.23k
  return sol;
4354
0
error:
4355
0
  sol_free(sol);
4356
0
  isl_basic_map_free(bmap);
4357
0
  return NULL;
4358
7.23k
}
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.48k
{
4369
3.48k
  isl_map *result = NULL;
4370
3.48k
  struct isl_sol *sol;
4371
3.48k
  struct isl_sol_map *sol_map;
4372
3.48k
4373
3.48k
  sol = basic_map_partial_lexopt_base_sol(bmap, dom, empty, max,
4374
3.48k
            &sol_map_init);
4375
3.48k
  if (!sol)
4376
0
    return NULL;
4377
3.48k
  sol_map = (struct isl_sol_map *) sol;
4378
3.48k
4379
3.48k
  result = isl_map_copy(sol_map->map);
4380
3.48k
  if (empty)
4381
2.63k
    *empty = isl_set_copy(sol_map->empty);
4382
3.48k
  sol_free(&sol_map->sol);
4383
3.48k
  return result;
4384
3.48k
}
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.45k
{
4392
7.45k
  int i, j;
4393
7.45k
  isl_ctx *ctx;
4394
7.45k
  int *occurrences;
4395
7.45k
4396
7.45k
  if (!bmap)
4397
0
    return NULL;
4398
7.45k
  ctx = isl_basic_map_get_ctx(bmap);
4399
7.45k
  occurrences = isl_calloc_array(ctx, int, n);
4400
7.45k
  if (!occurrences)
4401
0
    return NULL;
4402
7.45k
4403
26.2k
  
for (i = 0; 7.45k
i < bmap->n_ineq;
++i18.7k
) {
4404
99.6k
    for (j = 0; j < n; 
++j80.8k
) {
4405
80.8k
      if (!isl_int_is_zero(bmap->ineq[i][1 + j]))
4406
80.8k
        
occurrences[j]++15.0k
;
4407
80.8k
    }
4408
18.7k
  }
4409
7.45k
4410
7.45k
  return occurrences;
4411
7.45k
}
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
8.42k
{
4420
8.42k
  int i;
4421
8.42k
4422
31.6k
  for (i = 0; i < n; 
++i23.1k
) {
4423
25.3k
    if (isl_int_is_zero(c[i]))
4424
25.3k
      
continue22.3k
;
4425
3.04k
    if (occurrences[i] != 1)
4426
2.20k
      return 0;
4427
25.3k
  }
4428
8.42k
4429
8.42k
  
return 16.22k
;
4430
8.42k
}
4431
4432
/* Do all of the "n" initial variables that occur in inequality constraint
4433
 * "ineq" of "bmap" only occur in that constraint?
4434
 */
4435
static int all_single_occurrence(__isl_keep isl_basic_map *bmap, int ineq,
4436
  int n)
4437
0
{
4438
0
  int i, j;
4439
0
4440
0
  for (i = 0; i < n; ++i) {
4441
0
    if (isl_int_is_zero(bmap->ineq[ineq][1 + i]))
4442
0
      continue;
4443
0
    for (j = 0; j < bmap->n_ineq; ++j) {
4444
0
      if (j == ineq)
4445
0
        continue;
4446
0
      if (!isl_int_is_zero(bmap->ineq[j][1 + i]))
4447
0
        return 0;
4448
0
    }
4449
0
  }
4450
0
4451
0
  return 1;
4452
0
}
4453
4454
/* Structure used during detection of parallel constraints.
4455
 * n_in: number of "input" variables: isl_dim_param + isl_dim_in
4456
 * n_out: number of "output" variables: isl_dim_out + isl_dim_div
4457
 * val: the coefficients of the output variables
4458
 */
4459
struct isl_constraint_equal_info {
4460
  unsigned n_in;
4461
  unsigned n_out;
4462
  isl_int *val;
4463
};
4464
4465
/* Check whether the coefficients of the output variables
4466
 * of the constraint in "entry" are equal to info->val.
4467
 */
4468
static int constraint_equal(const void *entry, const void *val)
4469
224
{
4470
224
  isl_int **row = (isl_int **)entry;
4471
224
  const struct isl_constraint_equal_info *info = val;
4472
224
4473
224
  return isl_seq_eq((*row) + 1 + info->n_in, info->val, info->n_out);
4474
224
}
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.45k
{
4489
7.45k
  int i;
4490
7.45k
  isl_ctx *ctx;
4491
7.45k
  int *occurrences = NULL;
4492
7.45k
  struct isl_hash_table *table = NULL;
4493
7.45k
  struct isl_hash_table_entry *entry;
4494
7.45k
  struct isl_constraint_equal_info info;
4495
7.45k
  unsigned n_out;
4496
7.45k
  unsigned n_div;
4497
7.45k
4498
7.45k
  ctx = isl_basic_map_get_ctx(bmap);
4499
7.45k
  table = isl_hash_table_alloc(ctx, bmap->n_ineq);
4500
7.45k
  if (!table)
4501
0
    goto error;
4502
7.45k
4503
7.45k
  info.n_in = isl_basic_map_dim(bmap, isl_dim_param) +
4504
7.45k
        isl_basic_map_dim(bmap, isl_dim_in);
4505
7.45k
  occurrences = count_occurrences(bmap, info.n_in);
4506
7.45k
  if (info.n_in && 
!occurrences6.60k
)
4507
0
    goto error;
4508
7.45k
  n_out = isl_basic_map_dim(bmap, isl_dim_out);
4509
7.45k
  n_div = isl_basic_map_dim(bmap, isl_dim_div);
4510
7.45k
  info.n_out = n_out + n_div;
4511
25.7k
  for (i = 0; i < bmap->n_ineq; 
++i18.3k
) {
4512
18.5k
    uint32_t hash;
4513
18.5k
4514
18.5k
    info.val = bmap->ineq[i] + 1 + info.n_in;
4515
18.5k
    if (isl_seq_first_non_zero(info.val, n_out) < 0)
4516
10.0k
      continue;
4517
8.52k
    if (isl_seq_first_non_zero(info.val + n_out, n_div) >= 0)
4518
100
      continue;
4519
8.42k
    if (!single_occurrence(info.n_in, bmap->ineq[i] + 1,
4520
8.42k
          occurrences))
4521
2.20k
      continue;
4522
6.22k
    hash = isl_seq_get_hash(info.val, info.n_out);
4523
6.22k
    entry = isl_hash_table_find(ctx, table, hash,
4524
6.22k
              constraint_equal, &info, 1);
4525
6.22k
    if (!entry)
4526
0
      goto error;
4527
6.22k
    if (entry->data)
4528
224
      break;
4529
6.00k
    entry->data = &bmap->ineq[i];
4530
6.00k
  }
4531
7.45k
4532
7.45k
  if (i < bmap->n_ineq) {
4533
224
    *first = ((isl_int **)entry->data) - bmap->ineq; 
4534
224
    *second = i;
4535
224
  }
4536
7.45k
4537
7.45k
  isl_hash_table_free(ctx, table);
4538
7.45k
  free(occurrences);
4539
7.45k
4540
7.45k
  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.45k
}
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
788
{
4558
788
  isl_ctx *ctx;
4559
788
  int j, k;
4560
788
4561
788
  ctx = isl_mat_get_ctx(var);
4562
788
4563
2.36k
  for (j = 0; j < var->n_row; 
++j1.57k
) {
4564
1.57k
    if (j == i)
4565
788
      continue;
4566
788
    k = isl_basic_set_alloc_inequality(bset);
4567
788
    if (k < 0)
4568
0
      goto error;
4569
788
    isl_seq_combine(bset->ineq[k], ctx->one, var->row[j],
4570
788
        ctx->negone, var->row[i], var->n_col);
4571
788
    isl_int_set_si(bset->ineq[k][var->n_col], 0);
4572
788
    if (j < i)
4573
788
      
isl_int_sub_ui394
(bset->ineq[k][0], bset->ineq[k][0], 1);
4574
1.57k
  }
4575
788
4576
788
  bset = isl_basic_set_finalize(bset);
4577
788
4578
788
  return bset;
4579
0
error:
4580
0
  isl_basic_set_free(bset);
4581
0
  return NULL;
4582
788
}
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
224
{
4600
224
  int i, k;
4601
224
  isl_basic_set *bset = NULL;
4602
224
  isl_set *set = NULL;
4603
224
4604
224
  if (!dim || !var)
4605
0
    goto error;
4606
224
4607
224
  set = isl_set_alloc_space(isl_space_copy(dim),
4608
224
        var->n_row, ISL_SET_DISJOINT);
4609
224
4610
672
  for (i = 0; i < var->n_row; 
++i448
) {
4611
448
    bset = isl_basic_set_alloc_space(isl_space_copy(dim), 0,
4612
448
                 1, var->n_row - 1);
4613
448
    k = isl_basic_set_alloc_equality(bset);
4614
448
    if (k < 0)
4615
0
      goto error;
4616
448
    isl_seq_cpy(bset->eq[k], var->row[i], var->n_col);
4617
448
    isl_int_set_si(bset->eq[k][var->n_col], -1);
4618
448
    bset = select_minimum(bset, var, i);
4619
448
    set = isl_set_add_basic_set(set, bset);
4620
448
  }
4621
224
4622
224
  isl_space_free(dim);
4623
224
  isl_mat_free(var);
4624
224
  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
224
}
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
115
{
4645
115
  int i, j;
4646
115
  unsigned total;
4647
115
  unsigned pos;
4648
115
4649
115
  pos = cst->n_col - 1;
4650
115
  total = isl_basic_map_dim(bmap, isl_dim_all);
4651
115
4652
150
  for (i = 0; i < bmap->n_div; 
++i35
)
4653
115
    
if (113
!113
isl_int_is_zero113
(bmap->div[i][2 + pos]))
4654
113
      
return isl_bool_true78
;
4655
115
4656
115
  
for (i = 0; 37
i < bmap->n_eq76
;
++i39
)
4657
43
    if (!isl_int_is_zero(bmap->eq[i][1 + pos]))
4658
43
      
return isl_bool_true4
;
4659
37
4660
130
  
for (i = 0; 33
i < bmap->n_ineq;
++i97
) {
4661
118
    if (isl_int_is_nonneg(bmap->ineq[i][1 + pos]))
4662
118
      
continue61
;
4663
57
    if (!isl_int_is_negone(bmap->ineq[i][1 + pos]))
4664
57
      
return isl_bool_true0
;
4665
57
    if (isl_seq_first_non_zero(bmap->ineq[i] + 1 + pos + 1,
4666
57
             total - pos - 1) >= 0)
4667
15
      return isl_bool_true;
4668
42
4669
72
    
for (j = 0; 42
j < cst->n_row;
++j30
)
4670
66
      if (isl_seq_eq(bmap->ineq[i], cst->row[j], cst->n_col))
4671
36
        break;
4672
42
    if (j >= cst->n_row)
4673
6
      return isl_bool_true;
4674
118
  }
4675
33
4676
33
  
return isl_bool_false12
;
4677
115
}
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
12
{
4690
12
  return need_split_basic_map(bset_to_bmap(bset), cst);
4691
12
}
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
6
{
4699
6
  int i;
4700
6
4701
12
  for (i = 0; i < set->n; 
++i6
) {
4702
6
    isl_bool split;
4703
6
4704
6
    split = need_split_basic_set(set->p[i], cst);
4705
6
    if (split < 0 || split)
4706
0
      return split;
4707
6
  }
4708
6
4709
6
  return isl_bool_false;
4710
6
}
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 || !min_expr || !cst)
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->n; 
++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 (split)
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
54
{
4778
54
  int n_in;
4779
54
  int i;
4780
54
  isl_space *dim;
4781
54
  isl_map *res;
4782
54
4783
54
  if (!opt || !min_expr || !cst)
4784
0
    goto error;
4785
54
4786
54
  n_in = isl_map_dim(opt, isl_dim_in);
4787
54
  dim = isl_map_get_space(opt);
4788
54
  dim = isl_space_drop_dims(dim, isl_dim_in, n_in - 1, 1);
4789
54
  res = isl_map_empty(dim);
4790
54
4791
157
  for (i = 0; i < opt->n; 
++i103
) {
4792
103
    isl_map *map;
4793
103
    isl_bool split;
4794
103
4795
103
    map = isl_map_from_basic_map(isl_basic_map_copy(opt->p[i]));
4796
103
    split = need_split_basic_map(opt->p[i], cst);
4797
103
    if (split < 0)
4798
0
      map = isl_map_free(map);
4799
103
    else if (split)
4800
97
      map = isl_map_intersect_domain(map,
4801
97
                   isl_set_copy(min_expr));
4802
103
    map = isl_map_remove_dims(map, isl_dim_in, n_in - 1, 1);
4803
103
4804
103
    res = isl_map_union_disjoint(res, map);
4805
103
  }
4806
54
4807
54
  isl_map_free(opt);
4808
54
  isl_set_free(min_expr);
4809
54
  isl_mat_free(cst);
4810
54
  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
54
}
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
54
{
4836
54
  isl_map *opt;
4837
54
  isl_set *min_expr;
4838
54
4839
54
  min_expr = set_minimum(isl_basic_set_get_space(dom), isl_mat_copy(cst));
4840
54
4841
54
  opt = basic_map_partial_lexopt(bmap, dom, empty, max);
4842
54
4843
54
  if (empty) {
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
54
4849
54
  opt = split_domain(opt, min_expr, cst);
4850
54
  opt = isl_map_reset_space(opt, map_space);
4851
54
4852
54
  return opt;
4853
54
}
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.48k
{
4880
3.48k
  int n_div;
4881
3.48k
  int n_out;
4882
3.48k
4883
3.48k
  n_div = isl_basic_map_dim(bmap, isl_dim_div);
4884
3.48k
  n_out = isl_basic_map_dim(bmap, isl_dim_out);
4885
3.48k
  bmap = isl_basic_map_copy(bmap);
4886
3.48k
  if (ISL_FL_ISSET(flags, ISL_OPT_QE)) {
4887
180
    bmap = isl_basic_map_drop_constraints_involving_dims(bmap,
4888
180
              isl_dim_div, 0, n_div);
4889
180
    bmap = isl_basic_map_drop_constraints_involving_dims(bmap,
4890
180
              isl_dim_out, 0, n_out);
4891
180
  }
4892
3.48k
  return isl_basic_map_domain(bmap);
4893
3.48k
}
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
753
{
4910
753
}
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
754
{
4928
754
  int i, n;
4929
754
  isl_ctx *ctx;
4930
754
  isl_aff *aff;
4931
754
  isl_aff_list *list;
4932
754
4933
754
  if (sol->sol.error || !dom || !ma)
4934
0
    goto error;
4935
754
4936
754
  ctx = isl_basic_set_get_ctx(dom);
4937
754
  n = isl_multi_aff_dim(ma, isl_dim_out);
4938
754
  list = isl_aff_list_alloc(ctx, n);
4939
1.50k
  for (i = 0; i < n; 
++i754
) {
4940
754
    aff = isl_multi_aff_get_aff(ma, i);
4941
754
    list = isl_aff_list_add(list, aff);
4942
754
  }
4943
754
4944
754
  dom = isl_basic_set_finalize(dom);
4945
754
4946
754
  if (sol->fn(isl_basic_set_copy(dom), list, sol->user) < 0)
4947
0
    goto error;
4948
754
4949
754
  isl_basic_set_free(dom);
4950
754
  isl_multi_aff_free(ma);
4951
754
  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
754
{
4961
754
  sol_for_add((struct isl_sol_for *)sol, dom, ma);
4962
754
}
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
753
{
4969
753
  struct isl_sol_for *sol_for = NULL;
4970
753
  isl_space *dom_dim;
4971
753
  struct isl_basic_set *dom = NULL;
4972
753
4973
753
  sol_for = isl_calloc_type(bmap->ctx, struct isl_sol_for);
4974
753
  if (!sol_for)
4975
0
    goto error;
4976
753
4977
753
  dom_dim = isl_space_domain(isl_space_copy(bmap->dim));
4978
753
  dom = isl_basic_set_universe(dom_dim);
4979
753
4980
753
  sol_for->sol.free = &sol_for_free;
4981
753
  if (sol_init(&sol_for->sol, bmap, dom, max) < 0)
4982
0
    goto error;
4983
753
  sol_for->fn = fn;
4984
753
  sol_for->user = user;
4985
753
  sol_for->sol.add = &sol_for_add_wrap;
4986
753
  sol_for->sol.add_empty = NULL;
4987
753
4988
753
  isl_basic_set_free(dom);
4989
753
  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