MLIR  16.0.0git
Simplex.cpp
Go to the documentation of this file.
1 //===- Simplex.cpp - MLIR Simplex Class -----------------------------------===//
2 //
3 // Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
4 // See https://llvm.org/LICENSE.txt for license information.
5 // SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
6 //
7 //===----------------------------------------------------------------------===//
8 
12 #include "llvm/ADT/Optional.h"
13 #include "llvm/Support/Compiler.h"
14 #include <numeric>
15 
16 using namespace mlir;
17 using namespace presburger;
18 
20 
22 
23 // Return a + scale*b;
24 LLVM_ATTRIBUTE_UNUSED
27  assert(a.size() == b.size());
29  res.reserve(a.size());
30  for (unsigned i = 0, e = a.size(); i < e; ++i)
31  res.push_back(a[i] + scale * b[i]);
32  return res;
33 }
34 
35 SimplexBase::SimplexBase(unsigned nVar, bool mustUseBigM)
36  : usingBigM(mustUseBigM), nRedundant(0), nSymbol(0),
37  tableau(0, getNumFixedCols() + nVar), empty(false) {
38  colUnknown.insert(colUnknown.begin(), getNumFixedCols(), nullIndex);
39  for (unsigned i = 0; i < nVar; ++i) {
40  var.emplace_back(Orientation::Column, /*restricted=*/false,
41  /*pos=*/getNumFixedCols() + i);
42  colUnknown.push_back(i);
43  }
44 }
45 
46 SimplexBase::SimplexBase(unsigned nVar, bool mustUseBigM,
47  const llvm::SmallBitVector &isSymbol)
48  : SimplexBase(nVar, mustUseBigM) {
49  assert(isSymbol.size() == nVar && "invalid bitmask!");
50  // Invariant: nSymbol is the number of symbols that have been marked
51  // already and these occupy the columns
52  // [getNumFixedCols(), getNumFixedCols() + nSymbol).
53  for (unsigned symbolIdx : isSymbol.set_bits()) {
54  var[symbolIdx].isSymbol = true;
55  swapColumns(var[symbolIdx].pos, getNumFixedCols() + nSymbol);
56  ++nSymbol;
57  }
58 }
59 
61  assert(index != nullIndex && "nullIndex passed to unknownFromIndex");
62  return index >= 0 ? var[index] : con[~index];
63 }
64 
66  assert(col < getNumColumns() && "Invalid column");
67  return unknownFromIndex(colUnknown[col]);
68 }
69 
70 const Simplex::Unknown &SimplexBase::unknownFromRow(unsigned row) const {
71  assert(row < getNumRows() && "Invalid row");
72  return unknownFromIndex(rowUnknown[row]);
73 }
74 
76  assert(index != nullIndex && "nullIndex passed to unknownFromIndex");
77  return index >= 0 ? var[index] : con[~index];
78 }
79 
81  assert(col < getNumColumns() && "Invalid column");
82  return unknownFromIndex(colUnknown[col]);
83 }
84 
86  assert(row < getNumRows() && "Invalid row");
87  return unknownFromIndex(rowUnknown[row]);
88 }
89 
90 unsigned SimplexBase::addZeroRow(bool makeRestricted) {
91  // Resize the tableau to accommodate the extra row.
92  unsigned newRow = tableau.appendExtraRow();
93  assert(getNumRows() == getNumRows() && "Inconsistent tableau size");
94  rowUnknown.push_back(~con.size());
95  con.emplace_back(Orientation::Row, makeRestricted, newRow);
97  tableau(newRow, 0) = 1;
98  return newRow;
99 }
100 
101 /// Add a new row to the tableau corresponding to the given constant term and
102 /// list of coefficients. The coefficients are specified as a vector of
103 /// (variable index, coefficient) pairs.
104 unsigned SimplexBase::addRow(ArrayRef<MPInt> coeffs, bool makeRestricted) {
105  assert(coeffs.size() == var.size() + 1 &&
106  "Incorrect number of coefficients!");
107  assert(var.size() + getNumFixedCols() == getNumColumns() &&
108  "inconsistent column count!");
109 
110  unsigned newRow = addZeroRow(makeRestricted);
111  tableau(newRow, 1) = coeffs.back();
112  if (usingBigM) {
113  // When the lexicographic pivot rule is used, instead of the variables
114  //
115  // x, y, z ...
116  //
117  // we internally use the variables
118  //
119  // M, M + x, M + y, M + z, ...
120  //
121  // where M is the big M parameter. As such, when the user tries to add
122  // a row ax + by + cz + d, we express it in terms of our internal variables
123  // as -(a + b + c)M + a(M + x) + b(M + y) + c(M + z) + d.
124  //
125  // Symbols don't use the big M parameter since they do not get lex
126  // optimized.
127  MPInt bigMCoeff(0);
128  for (unsigned i = 0; i < coeffs.size() - 1; ++i)
129  if (!var[i].isSymbol)
130  bigMCoeff -= coeffs[i];
131  // The coefficient to the big M parameter is stored in column 2.
132  tableau(newRow, 2) = bigMCoeff;
133  }
134 
135  // Process each given variable coefficient.
136  for (unsigned i = 0; i < var.size(); ++i) {
137  unsigned pos = var[i].pos;
138  if (coeffs[i] == 0)
139  continue;
140 
141  if (var[i].orientation == Orientation::Column) {
142  // If a variable is in column position at column col, then we just add the
143  // coefficient for that variable (scaled by the common row denominator) to
144  // the corresponding entry in the new row.
145  tableau(newRow, pos) += coeffs[i] * tableau(newRow, 0);
146  continue;
147  }
148 
149  // If the variable is in row position, we need to add that row to the new
150  // row, scaled by the coefficient for the variable, accounting for the two
151  // rows potentially having different denominators. The new denominator is
152  // the lcm of the two.
153  MPInt lcm = presburger::lcm(tableau(newRow, 0), tableau(pos, 0));
154  MPInt nRowCoeff = lcm / tableau(newRow, 0);
155  MPInt idxRowCoeff = coeffs[i] * (lcm / tableau(pos, 0));
156  tableau(newRow, 0) = lcm;
157  for (unsigned col = 1, e = getNumColumns(); col < e; ++col)
158  tableau(newRow, col) =
159  nRowCoeff * tableau(newRow, col) + idxRowCoeff * tableau(pos, col);
160  }
161 
162  tableau.normalizeRow(newRow);
163  // Push to undo log along with the index of the new constraint.
164  return con.size() - 1;
165 }
166 
167 namespace {
168 bool signMatchesDirection(const MPInt &elem, Direction direction) {
169  assert(elem != 0 && "elem should not be 0");
170  return direction == Direction::Up ? elem > 0 : elem < 0;
171 }
172 
173 Direction flippedDirection(Direction direction) {
174  return direction == Direction::Up ? Direction::Down : Simplex::Direction::Up;
175 }
176 } // namespace
177 
178 /// We simply make the tableau consistent while maintaining a lexicopositive
179 /// basis transform, and then return the sample value. If the tableau becomes
180 /// empty, we return empty.
181 ///
182 /// Let the variables be x = (x_1, ... x_n).
183 /// Let the basis unknowns be y = (y_1, ... y_n).
184 /// We have that x = A*y + b for some n x n matrix A and n x 1 column vector b.
185 ///
186 /// As we will show below, A*y is either zero or lexicopositive.
187 /// Adding a lexicopositive vector to b will make it lexicographically
188 /// greater, so A*y + b is always equal to or lexicographically greater than b.
189 /// Thus, since we can attain x = b, that is the lexicographic minimum.
190 ///
191 /// We have that that every column in A is lexicopositive, i.e., has at least
192 /// one non-zero element, with the first such element being positive. Since for
193 /// the tableau to be consistent we must have non-negative sample values not
194 /// only for the constraints but also for the variables, we also have x >= 0 and
195 /// y >= 0, by which we mean every element in these vectors is non-negative.
196 ///
197 /// Proof that if every column in A is lexicopositive, and y >= 0, then
198 /// A*y is zero or lexicopositive. Begin by considering A_1, the first row of A.
199 /// If this row is all zeros, then (A*y)_1 = (A_1)*y = 0; proceed to the next
200 /// row. If we run out of rows, A*y is zero and we are done; otherwise, we
201 /// encounter some row A_i that has a non-zero element. Every column is
202 /// lexicopositive and so has some positive element before any negative elements
203 /// occur, so the element in this row for any column, if non-zero, must be
204 /// positive. Consider (A*y)_i = (A_i)*y. All the elements in both vectors are
205 /// non-negative, so if this is non-zero then it must be positive. Then the
206 /// first non-zero element of A*y is positive so A*y is lexicopositive.
207 ///
208 /// Otherwise, if (A_i)*y is zero, then for every column j that had a non-zero
209 /// element in A_i, y_j is zero. Thus these columns have no contribution to A*y
210 /// and we can completely ignore these columns of A. We now continue downwards,
211 /// looking for rows of A that have a non-zero element other than in the ignored
212 /// columns. If we find one, say A_k, once again these elements must be positive
213 /// since they are the first non-zero element in each of these columns, so if
214 /// (A_k)*y is not zero then we have that A*y is lexicopositive and if not we
215 /// add these to the set of ignored columns and continue to the next row. If we
216 /// run out of rows, then A*y is zero and we are done.
218  if (restoreRationalConsistency().failed()) {
219  markEmpty();
220  return OptimumKind::Empty;
221  }
222  return getRationalSample();
223 }
224 
225 /// Given a row that has a non-integer sample value, add an inequality such
226 /// that this fractional sample value is cut away from the polytope. The added
227 /// inequality will be such that no integer points are removed. i.e., the
228 /// integer lexmin, if it exists, is the same with and without this constraint.
229 ///
230 /// Let the row be
231 /// (c + coeffM*M + a_1*s_1 + ... + a_m*s_m + b_1*y_1 + ... + b_n*y_n)/d,
232 /// where s_1, ... s_m are the symbols and
233 /// y_1, ... y_n are the other basis unknowns.
234 ///
235 /// For this to be an integer, we want
236 /// coeffM*M + a_1*s_1 + ... + a_m*s_m + b_1*y_1 + ... + b_n*y_n = -c (mod d)
237 /// Note that this constraint must always hold, independent of the basis,
238 /// becuse the row unknown's value always equals this expression, even if *we*
239 /// later compute the sample value from a different expression based on a
240 /// different basis.
241 ///
242 /// Let us assume that M has a factor of d in it. Imposing this constraint on M
243 /// does not in any way hinder us from finding a value of M that is big enough.
244 /// Moreover, this function is only called when the symbolic part of the sample,
245 /// a_1*s_1 + ... + a_m*s_m, is known to be an integer.
246 ///
247 /// Also, we can safely reduce the coefficients modulo d, so we have:
248 ///
249 /// (b_1%d)y_1 + ... + (b_n%d)y_n = (-c%d) + k*d for some integer `k`
250 ///
251 /// Note that all coefficient modulos here are non-negative. Also, all the
252 /// unknowns are non-negative here as both constraints and variables are
253 /// non-negative in LexSimplexBase. (We used the big M trick to make the
254 /// variables non-negative). Therefore, the LHS here is non-negative.
255 /// Since 0 <= (-c%d) < d, k is the quotient of dividing the LHS by d and
256 /// is therefore non-negative as well.
257 ///
258 /// So we have
259 /// ((b_1%d)y_1 + ... + (b_n%d)y_n - (-c%d))/d >= 0.
260 ///
261 /// The constraint is violated when added (it would be useless otherwise)
262 /// so we immediately try to move it to a column.
264  MPInt d = tableau(row, 0);
265  unsigned cutRow = addZeroRow(/*makeRestricted=*/true);
266  tableau(cutRow, 0) = d;
267  tableau(cutRow, 1) = -mod(-tableau(row, 1), d); // -c%d.
268  tableau(cutRow, 2) = 0;
269  for (unsigned col = 3 + nSymbol, e = getNumColumns(); col < e; ++col)
270  tableau(cutRow, col) = mod(tableau(row, col), d); // b_i%d.
271  return moveRowUnknownToColumn(cutRow);
272 }
273 
274 Optional<unsigned> LexSimplex::maybeGetNonIntegralVarRow() const {
275  for (const Unknown &u : var) {
276  if (u.orientation == Orientation::Column)
277  continue;
278  // If the sample value is of the form (a/d)M + b/d, we need b to be
279  // divisible by d. We assume M contains all possible
280  // factors and is divisible by everything.
281  unsigned row = u.pos;
282  if (tableau(row, 1) % tableau(row, 0) != 0)
283  return row;
284  }
285  return {};
286 }
287 
289  // We first try to make the tableau consistent.
290  if (restoreRationalConsistency().failed())
291  return OptimumKind::Empty;
292 
293  // Then, if the sample value is integral, we are done.
294  while (Optional<unsigned> maybeRow = maybeGetNonIntegralVarRow()) {
295  // Otherwise, for the variable whose row has a non-integral sample value,
296  // we add a cut, a constraint that remove this rational point
297  // while preserving all integer points, thus keeping the lexmin the same.
298  // We then again try to make the tableau with the new constraint
299  // consistent. This continues until the tableau becomes empty, in which
300  // case there is no integer point, or until there are no variables with
301  // non-integral sample values.
302  //
303  // Failure indicates that the tableau became empty, which occurs when the
304  // polytope is integer empty.
305  if (addCut(*maybeRow).failed())
306  return OptimumKind::Empty;
307  if (restoreRationalConsistency().failed())
308  return OptimumKind::Empty;
309  }
310 
311  MaybeOptimum<SmallVector<Fraction, 8>> sample = getRationalSample();
312  assert(!sample.isEmpty() && "If we reached here the sample should exist!");
313  if (sample.isUnbounded())
314  return OptimumKind::Unbounded;
315  return llvm::to_vector<8>(
316  llvm::map_range(*sample, std::mem_fn(&Fraction::getAsInteger)));
317 }
318 
320  SimplexRollbackScopeExit scopeExit(*this);
321  addInequality(coeffs);
322  return findIntegerLexMin().isEmpty();
323 }
324 
326  return isSeparateInequality(getComplementIneq(coeffs));
327 }
328 
330 SymbolicLexSimplex::getSymbolicSampleNumerator(unsigned row) const {
331  SmallVector<MPInt, 8> sample;
332  sample.reserve(nSymbol + 1);
333  for (unsigned col = 3; col < 3 + nSymbol; ++col)
334  sample.push_back(tableau(row, col));
335  sample.push_back(tableau(row, 1));
336  return sample;
337 }
338 
340 SymbolicLexSimplex::getSymbolicSampleIneq(unsigned row) const {
341  SmallVector<MPInt, 8> sample = getSymbolicSampleNumerator(row);
342  // The inequality is equivalent to the GCD-normalized one.
343  normalizeRange(sample);
344  return sample;
345 }
346 
348  appendVariable();
349  swapColumns(3 + nSymbol, getNumColumns() - 1);
350  var.back().isSymbol = true;
351  nSymbol++;
352 }
353 
354 static bool isRangeDivisibleBy(ArrayRef<MPInt> range, const MPInt &divisor) {
355  assert(divisor > 0 && "divisor must be positive!");
356  return llvm::all_of(range,
357  [divisor](const MPInt &x) { return x % divisor == 0; });
358 }
359 
360 bool SymbolicLexSimplex::isSymbolicSampleIntegral(unsigned row) const {
361  MPInt denom = tableau(row, 0);
362  return tableau(row, 1) % denom == 0 &&
363  isRangeDivisibleBy(tableau.getRow(row).slice(3, nSymbol), denom);
364 }
365 
366 /// This proceeds similarly to LexSimplexBase::addCut(). We are given a row that
367 /// has a symbolic sample value with fractional coefficients.
368 ///
369 /// Let the row be
370 /// (c + coeffM*M + sum_i a_i*s_i + sum_j b_j*y_j)/d,
371 /// where s_1, ... s_m are the symbols and
372 /// y_1, ... y_n are the other basis unknowns.
373 ///
374 /// As in LexSimplex::addCut, for this to be an integer, we want
375 ///
376 /// coeffM*M + sum_j b_j*y_j = -c + sum_i (-a_i*s_i) (mod d)
377 ///
378 /// This time, a_1*s_1 + ... + a_m*s_m may not be an integer. We find that
379 ///
380 /// sum_i (b_i%d)y_i = ((-c%d) + sum_i (-a_i%d)s_i)%d + k*d for some integer k
381 ///
382 /// where we take a modulo of the whole symbolic expression on the right to
383 /// bring it into the range [0, d - 1]. Therefore, as in addCut(),
384 /// k is the quotient on dividing the LHS by d, and since LHS >= 0, we have
385 /// k >= 0 as well. If all the a_i are divisible by d, then we can add the
386 /// constraint directly. Otherwise, we realize the modulo of the symbolic
387 /// expression by adding a division variable
388 ///
389 /// q = ((-c%d) + sum_i (-a_i%d)s_i)/d
390 ///
391 /// to the symbol domain, so the equality becomes
392 ///
393 /// sum_i (b_i%d)y_i = (-c%d) + sum_i (-a_i%d)s_i - q*d + k*d for some integer k
394 ///
395 /// So the cut is
396 /// (sum_i (b_i%d)y_i - (-c%d) - sum_i (-a_i%d)s_i + q*d)/d >= 0
397 /// This constraint is violated when added so we immediately try to move it to a
398 /// column.
399 LogicalResult SymbolicLexSimplex::addSymbolicCut(unsigned row) {
400  MPInt d = tableau(row, 0);
401  if (isRangeDivisibleBy(tableau.getRow(row).slice(3, nSymbol), d)) {
402  // The coefficients of symbols in the symbol numerator are divisible
403  // by the denominator, so we can add the constraint directly,
404  // i.e., ignore the symbols and add a regular cut as in addCut().
405  return addCut(row);
406  }
407 
408  // Construct the division variable `q = ((-c%d) + sum_i (-a_i%d)s_i)/d`.
409  SmallVector<MPInt, 8> divCoeffs;
410  divCoeffs.reserve(nSymbol + 1);
411  MPInt divDenom = d;
412  for (unsigned col = 3; col < 3 + nSymbol; ++col)
413  divCoeffs.push_back(mod(-tableau(row, col), divDenom)); // (-a_i%d)s_i
414  divCoeffs.push_back(mod(-tableau(row, 1), divDenom)); // -c%d.
415  normalizeDiv(divCoeffs, divDenom);
416 
417  domainSimplex.addDivisionVariable(divCoeffs, divDenom);
418  domainPoly.addLocalFloorDiv(divCoeffs, divDenom);
419 
420  // Update `this` to account for the additional symbol we just added.
421  appendSymbol();
422 
423  // Add the cut (sum_i (b_i%d)y_i - (-c%d) + sum_i -(-a_i%d)s_i + q*d)/d >= 0.
424  unsigned cutRow = addZeroRow(/*makeRestricted=*/true);
425  tableau(cutRow, 0) = d;
426  tableau(cutRow, 2) = 0;
427 
428  tableau(cutRow, 1) = -mod(-tableau(row, 1), d); // -(-c%d).
429  for (unsigned col = 3; col < 3 + nSymbol - 1; ++col)
430  tableau(cutRow, col) = -mod(-tableau(row, col), d); // -(-a_i%d)s_i.
431  tableau(cutRow, 3 + nSymbol - 1) = d; // q*d.
432 
433  for (unsigned col = 3 + nSymbol, e = getNumColumns(); col < e; ++col)
434  tableau(cutRow, col) = mod(tableau(row, col), d); // (b_i%d)y_i.
435  return moveRowUnknownToColumn(cutRow);
436 }
437 
438 void SymbolicLexSimplex::recordOutput(SymbolicLexMin &result) const {
439  Matrix output(0, domainPoly.getNumVars() + 1);
440  output.reserveRows(result.lexmin.getNumOutputs());
441  for (const Unknown &u : var) {
442  if (u.isSymbol)
443  continue;
444 
445  if (u.orientation == Orientation::Column) {
446  // M + u has a sample value of zero so u has a sample value of -M, i.e,
447  // unbounded.
448  result.unboundedDomain.unionInPlace(domainPoly);
449  return;
450  }
451 
452  MPInt denom = tableau(u.pos, 0);
453  if (tableau(u.pos, 2) < denom) {
454  // M + u has a sample value of fM + something, where f < 1, so
455  // u = (f - 1)M + something, which has a negative coefficient for M,
456  // and so is unbounded.
457  result.unboundedDomain.unionInPlace(domainPoly);
458  return;
459  }
460  assert(tableau(u.pos, 2) == denom &&
461  "Coefficient of M should not be greater than 1!");
462 
463  SmallVector<MPInt, 8> sample = getSymbolicSampleNumerator(u.pos);
464  for (MPInt &elem : sample) {
465  assert(elem % denom == 0 && "coefficients must be integral!");
466  elem /= denom;
467  }
468  output.appendExtraRow(sample);
469  }
470 
471  // Store the output in a MultiAffineFunction and add it the result.
472  PresburgerSpace funcSpace = result.lexmin.getSpace();
473  funcSpace.insertVar(VarKind::Local, 0, domainPoly.getNumLocalVars());
474 
475  result.lexmin.addPiece(
476  {PresburgerSet(domainPoly),
477  MultiAffineFunction(funcSpace, output, domainPoly.getLocalReprs())});
478 }
479 
480 Optional<unsigned> SymbolicLexSimplex::maybeGetAlwaysViolatedRow() {
481  // First look for rows that are clearly violated just from the big M
482  // coefficient, without needing to perform any simplex queries on the domain.
483  for (unsigned row = 0, e = getNumRows(); row < e; ++row)
484  if (tableau(row, 2) < 0)
485  return row;
486 
487  for (unsigned row = 0, e = getNumRows(); row < e; ++row) {
488  if (tableau(row, 2) > 0)
489  continue;
490  if (domainSimplex.isSeparateInequality(getSymbolicSampleIneq(row))) {
491  // Sample numerator always takes negative values in the symbol domain.
492  return row;
493  }
494  }
495  return {};
496 }
497 
498 Optional<unsigned> SymbolicLexSimplex::maybeGetNonIntegralVarRow() {
499  for (const Unknown &u : var) {
500  if (u.orientation == Orientation::Column)
501  continue;
502  assert(!u.isSymbol && "Symbol should not be in row orientation!");
503  if (!isSymbolicSampleIntegral(u.pos))
504  return u.pos;
505  }
506  return {};
507 }
508 
509 /// The non-branching pivots are just the ones moving the rows
510 /// that are always violated in the symbol domain.
511 LogicalResult SymbolicLexSimplex::doNonBranchingPivots() {
512  while (Optional<unsigned> row = maybeGetAlwaysViolatedRow())
513  if (moveRowUnknownToColumn(*row).failed())
514  return failure();
515  return success();
516 }
517 
520  /*numDomain=*/domainPoly.getNumDimVars(),
521  /*numRange=*/var.size() - nSymbol,
522  /*numSymbols=*/domainPoly.getNumSymbolVars()));
523 
524  /// The algorithm is more naturally expressed recursively, but we implement
525  /// it iteratively here to avoid potential issues with stack overflows in the
526  /// compiler. We explicitly maintain the stack frames in a vector.
527  ///
528  /// To "recurse", we store the current "stack frame", i.e., state variables
529  /// that we will need when we "return", into `stack`, increment `level`, and
530  /// `continue`. To "tail recurse", we just `continue`.
531  /// To "return", we decrement `level` and `continue`.
532  ///
533  /// When there is no stack frame for the current `level`, this indicates that
534  /// we have just "recursed" or "tail recursed". When there does exist one,
535  /// this indicates that we have just "returned" from recursing. There is only
536  /// one point at which non-tail calls occur so we always "return" there.
537  unsigned level = 1;
538  struct StackFrame {
539  int splitIndex;
540  unsigned snapshot;
541  unsigned domainSnapshot;
542  IntegerRelation::CountsSnapshot domainPolyCounts;
543  };
545 
546  while (level > 0) {
547  assert(level >= stack.size());
548  if (level > stack.size()) {
549  if (empty || domainSimplex.findIntegerLexMin().isEmpty()) {
550  // No integer points; return.
551  --level;
552  continue;
553  }
554 
555  if (doNonBranchingPivots().failed()) {
556  // Could not find pivots for violated constraints; return.
557  --level;
558  continue;
559  }
560 
561  SmallVector<MPInt, 8> symbolicSample;
562  unsigned splitRow = 0;
563  for (unsigned e = getNumRows(); splitRow < e; ++splitRow) {
564  if (tableau(splitRow, 2) > 0)
565  continue;
566  assert(tableau(splitRow, 2) == 0 &&
567  "Non-branching pivots should have been handled already!");
568 
569  symbolicSample = getSymbolicSampleIneq(splitRow);
570  if (domainSimplex.isRedundantInequality(symbolicSample))
571  continue;
572 
573  // It's neither redundant nor separate, so it takes both positive and
574  // negative values, and hence constitutes a row for which we need to
575  // split the domain and separately run each case.
576  assert(!domainSimplex.isSeparateInequality(symbolicSample) &&
577  "Non-branching pivots should have been handled already!");
578  break;
579  }
580 
581  if (splitRow < getNumRows()) {
582  unsigned domainSnapshot = domainSimplex.getSnapshot();
583  IntegerRelation::CountsSnapshot domainPolyCounts =
584  domainPoly.getCounts();
585 
586  // First, we consider the part of the domain where the row is not
587  // violated. We don't have to do any pivots for the row in this case,
588  // but we record the additional constraint that defines this part of
589  // the domain.
590  domainSimplex.addInequality(symbolicSample);
591  domainPoly.addInequality(symbolicSample);
592 
593  // Recurse.
594  //
595  // On return, the basis as a set is preserved but not the internal
596  // ordering within rows or columns. Thus, we take note of the index of
597  // the Unknown that caused the split, which may be in a different
598  // row when we come back from recursing. We will need this to recurse
599  // on the other part of the split domain, where the row is violated.
600  //
601  // Note that we have to capture the index above and not a reference to
602  // the Unknown itself, since the array it lives in might get
603  // reallocated.
604  int splitIndex = rowUnknown[splitRow];
605  unsigned snapshot = getSnapshot();
606  stack.push_back(
607  {splitIndex, snapshot, domainSnapshot, domainPolyCounts});
608  ++level;
609  continue;
610  }
611 
612  // The tableau is rationally consistent for the current domain.
613  // Now we look for non-integral sample values and add cuts for them.
614  if (Optional<unsigned> row = maybeGetNonIntegralVarRow()) {
615  if (addSymbolicCut(*row).failed()) {
616  // No integral points; return.
617  --level;
618  continue;
619  }
620 
621  // Rerun this level with the added cut constraint (tail recurse).
622  continue;
623  }
624 
625  // Record output and return.
626  recordOutput(result);
627  --level;
628  continue;
629  }
630 
631  if (level == stack.size()) {
632  // We have "returned" from "recursing".
633  const StackFrame &frame = stack.back();
634  domainPoly.truncate(frame.domainPolyCounts);
635  domainSimplex.rollback(frame.domainSnapshot);
636  rollback(frame.snapshot);
637  const Unknown &u = unknownFromIndex(frame.splitIndex);
638 
639  // Drop the frame. We don't need it anymore.
640  stack.pop_back();
641 
642  // Now we consider the part of the domain where the unknown `splitIndex`
643  // was negative.
644  assert(u.orientation == Orientation::Row &&
645  "The split row should have been returned to row orientation!");
646  SmallVector<MPInt, 8> splitIneq =
647  getComplementIneq(getSymbolicSampleIneq(u.pos));
648  normalizeRange(splitIneq);
649  if (moveRowUnknownToColumn(u.pos).failed()) {
650  // The unknown can't be made non-negative; return.
651  --level;
652  continue;
653  }
654 
655  // The unknown can be made negative; recurse with the corresponding domain
656  // constraints.
657  domainSimplex.addInequality(splitIneq);
658  domainPoly.addInequality(splitIneq);
659 
660  // We are now taking care of the second half of the domain and we don't
661  // need to do anything else here after returning, so it's a tail recurse.
662  continue;
663  }
664  }
665 
666  return result;
667 }
668 
669 bool LexSimplex::rowIsViolated(unsigned row) const {
670  if (tableau(row, 2) < 0)
671  return true;
672  if (tableau(row, 2) == 0 && tableau(row, 1) < 0)
673  return true;
674  return false;
675 }
676 
677 Optional<unsigned> LexSimplex::maybeGetViolatedRow() const {
678  for (unsigned row = 0, e = getNumRows(); row < e; ++row)
679  if (rowIsViolated(row))
680  return row;
681  return {};
682 }
683 
684 /// We simply look for violated rows and keep trying to move them to column
685 /// orientation, which always succeeds unless the constraints have no solution
686 /// in which case we just give up and return.
687 LogicalResult LexSimplex::restoreRationalConsistency() {
688  if (empty)
689  return failure();
690  while (Optional<unsigned> maybeViolatedRow = maybeGetViolatedRow())
691  if (moveRowUnknownToColumn(*maybeViolatedRow).failed())
692  return failure();
693  return success();
694 }
695 
696 // Move the row unknown to column orientation while preserving lexicopositivity
697 // of the basis transform. The sample value of the row must be non-positive.
698 //
699 // We only consider pivots where the pivot element is positive. Suppose no such
700 // pivot exists, i.e., some violated row has no positive coefficient for any
701 // basis unknown. The row can be represented as (s + c_1*u_1 + ... + c_n*u_n)/d,
702 // where d is the denominator, s is the sample value and the c_i are the basis
703 // coefficients. If s != 0, then since any feasible assignment of the basis
704 // satisfies u_i >= 0 for all i, and we have s < 0 as well as c_i < 0 for all i,
705 // any feasible assignment would violate this row and therefore the constraints
706 // have no solution.
707 //
708 // We can preserve lexicopositivity by picking the pivot column with positive
709 // pivot element that makes the lexicographically smallest change to the sample
710 // point.
711 //
712 // Proof. Let
713 // x = (x_1, ... x_n) be the variables,
714 // z = (z_1, ... z_m) be the constraints,
715 // y = (y_1, ... y_n) be the current basis, and
716 // define w = (x_1, ... x_n, z_1, ... z_m) = B*y + s.
717 // B is basically the simplex tableau of our implementation except that instead
718 // of only describing the transform to get back the non-basis unknowns, it
719 // defines the values of all the unknowns in terms of the basis unknowns.
720 // Similarly, s is the column for the sample value.
721 //
722 // Our goal is to show that each column in B, restricted to the first n
723 // rows, is lexicopositive after the pivot if it is so before. This is
724 // equivalent to saying the columns in the whole matrix are lexicopositive;
725 // there must be some non-zero element in every column in the first n rows since
726 // the n variables cannot be spanned without using all the n basis unknowns.
727 //
728 // Consider a pivot where z_i replaces y_j in the basis. Recall the pivot
729 // transform for the tableau derived for SimplexBase::pivot:
730 //
731 // pivot col other col pivot col other col
732 // pivot row a b -> pivot row 1/a -b/a
733 // other row c d other row c/a d - bc/a
734 //
735 // Similarly, a pivot results in B changing to B' and c to c'; the difference
736 // between the tableau and these matrices B and B' is that there is no special
737 // case for the pivot row, since it continues to represent the same unknown. The
738 // same formula applies for all rows:
739 //
740 // B'.col(j) = B.col(j) / B(i,j)
741 // B'.col(k) = B.col(k) - B(i,k) * B.col(j) / B(i,j) for k != j
742 // and similarly, s' = s - s_i * B.col(j) / B(i,j).
743 //
744 // If s_i == 0, then the sample value remains unchanged. Otherwise, if s_i < 0,
745 // the change in sample value when pivoting with column a is lexicographically
746 // smaller than that when pivoting with column b iff B.col(a) / B(i, a) is
747 // lexicographically smaller than B.col(b) / B(i, b).
748 //
749 // Since B(i, j) > 0, column j remains lexicopositive.
750 //
751 // For the other columns, suppose C.col(k) is not lexicopositive.
752 // This means that for some p, for all t < p,
753 // C(t,k) = 0 => B(t,k) = B(t,j) * B(i,k) / B(i,j) and
754 // C(t,k) < 0 => B(p,k) < B(t,j) * B(i,k) / B(i,j),
755 // which is in contradiction to the fact that B.col(j) / B(i,j) must be
756 // lexicographically smaller than B.col(k) / B(i,k), since it lexicographically
757 // minimizes the change in sample value.
759  Optional<unsigned> maybeColumn;
760  for (unsigned col = 3 + nSymbol, e = getNumColumns(); col < e; ++col) {
761  if (tableau(row, col) <= 0)
762  continue;
763  maybeColumn =
764  !maybeColumn ? col : getLexMinPivotColumn(row, *maybeColumn, col);
765  }
766 
767  if (!maybeColumn)
768  return failure();
769 
770  pivot(row, *maybeColumn);
771  return success();
772 }
773 
774 unsigned LexSimplexBase::getLexMinPivotColumn(unsigned row, unsigned colA,
775  unsigned colB) const {
776  // First, let's consider the non-symbolic case.
777  // A pivot causes the following change. (in the diagram the matrix elements
778  // are shown as rationals and there is no common denominator used)
779  //
780  // pivot col big M col const col
781  // pivot row a p b
782  // other row c q d
783  // |
784  // v
785  //
786  // pivot col big M col const col
787  // pivot row 1/a -p/a -b/a
788  // other row c/a q - pc/a d - bc/a
789  //
790  // Let the sample value of the pivot row be s = pM + b before the pivot. Since
791  // the pivot row represents a violated constraint we know that s < 0.
792  //
793  // If the variable is a non-pivot column, its sample value is zero before and
794  // after the pivot.
795  //
796  // If the variable is the pivot column, then its sample value goes from 0 to
797  // (-p/a)M + (-b/a), i.e. 0 to -(pM + b)/a. Thus the change in the sample
798  // value is -s/a.
799  //
800  // If the variable is the pivot row, its sample value goes from s to 0, for a
801  // change of -s.
802  //
803  // If the variable is a non-pivot row, its sample value changes from
804  // qM + d to qM + d + (-pc/a)M + (-bc/a). Thus the change in sample value
805  // is -(pM + b)(c/a) = -sc/a.
806  //
807  // Thus the change in sample value is either 0, -s/a, -s, or -sc/a. Here -s is
808  // fixed for all calls to this function since the row and tableau are fixed.
809  // The callee just wants to compare the return values with the return value of
810  // other invocations of the same function. So the -s is common for all
811  // comparisons involved and can be ignored, since -s is strictly positive.
812  //
813  // Thus we take away this common factor and just return 0, 1/a, 1, or c/a as
814  // appropriate. This allows us to run the entire algorithm treating M
815  // symbolically, as the pivot to be performed does not depend on the value
816  // of M, so long as the sample value s is negative. Note that this is not
817  // because of any special feature of M; by the same argument, we ignore the
818  // symbols too. The caller ensure that the sample value s is negative for
819  // all possible values of the symbols.
820  auto getSampleChangeCoeffForVar = [this, row](unsigned col,
821  const Unknown &u) -> Fraction {
822  MPInt a = tableau(row, col);
823  if (u.orientation == Orientation::Column) {
824  // Pivot column case.
825  if (u.pos == col)
826  return {1, a};
827 
828  // Non-pivot column case.
829  return {0, 1};
830  }
831 
832  // Pivot row case.
833  if (u.pos == row)
834  return {1, 1};
835 
836  // Non-pivot row case.
837  MPInt c = tableau(u.pos, col);
838  return {c, a};
839  };
840 
841  for (const Unknown &u : var) {
842  Fraction changeA = getSampleChangeCoeffForVar(colA, u);
843  Fraction changeB = getSampleChangeCoeffForVar(colB, u);
844  if (changeA < changeB)
845  return colA;
846  if (changeA > changeB)
847  return colB;
848  }
849 
850  // If we reached here, both result in exactly the same changes, so it
851  // doesn't matter which we return.
852  return colA;
853 }
854 
855 /// Find a pivot to change the sample value of the row in the specified
856 /// direction. The returned pivot row will involve `row` if and only if the
857 /// unknown is unbounded in the specified direction.
858 ///
859 /// To increase (resp. decrease) the value of a row, we need to find a live
860 /// column with a non-zero coefficient. If the coefficient is positive, we need
861 /// to increase (decrease) the value of the column, and if the coefficient is
862 /// negative, we need to decrease (increase) the value of the column. Also,
863 /// we cannot decrease the sample value of restricted columns.
864 ///
865 /// If multiple columns are valid, we break ties by considering a lexicographic
866 /// ordering where we prefer unknowns with lower index.
867 Optional<SimplexBase::Pivot> Simplex::findPivot(int row,
868  Direction direction) const {
869  Optional<unsigned> col;
870  for (unsigned j = 2, e = getNumColumns(); j < e; ++j) {
871  MPInt elem = tableau(row, j);
872  if (elem == 0)
873  continue;
874 
875  if (unknownFromColumn(j).restricted &&
876  !signMatchesDirection(elem, direction))
877  continue;
878  if (!col || colUnknown[j] < colUnknown[*col])
879  col = j;
880  }
881 
882  if (!col)
883  return {};
884 
885  Direction newDirection =
886  tableau(row, *col) < 0 ? flippedDirection(direction) : direction;
887  Optional<unsigned> maybePivotRow = findPivotRow(row, newDirection, *col);
888  return Pivot{maybePivotRow.value_or(row), *col};
889 }
890 
891 /// Swap the associated unknowns for the row and the column.
892 ///
893 /// First we swap the index associated with the row and column. Then we update
894 /// the unknowns to reflect their new position and orientation.
895 void SimplexBase::swapRowWithCol(unsigned row, unsigned col) {
896  std::swap(rowUnknown[row], colUnknown[col]);
897  Unknown &uCol = unknownFromColumn(col);
898  Unknown &uRow = unknownFromRow(row);
901  uCol.pos = col;
902  uRow.pos = row;
903 }
904 
905 void SimplexBase::pivot(Pivot pair) { pivot(pair.row, pair.column); }
906 
907 /// Pivot pivotRow and pivotCol.
908 ///
909 /// Let R be the pivot row unknown and let C be the pivot col unknown.
910 /// Since initially R = a*C + sum b_i * X_i
911 /// (where the sum is over the other column's unknowns, x_i)
912 /// C = (R - (sum b_i * X_i))/a
913 ///
914 /// Let u be some other row unknown.
915 /// u = c*C + sum d_i * X_i
916 /// So u = c*(R - sum b_i * X_i)/a + sum d_i * X_i
917 ///
918 /// This results in the following transform:
919 /// pivot col other col pivot col other col
920 /// pivot row a b -> pivot row 1/a -b/a
921 /// other row c d other row c/a d - bc/a
922 ///
923 /// Taking into account the common denominators p and q:
924 ///
925 /// pivot col other col pivot col other col
926 /// pivot row a/p b/p -> pivot row p/a -b/a
927 /// other row c/q d/q other row cp/aq (da - bc)/aq
928 ///
929 /// The pivot row transform is accomplished be swapping a with the pivot row's
930 /// common denominator and negating the pivot row except for the pivot column
931 /// element.
932 void SimplexBase::pivot(unsigned pivotRow, unsigned pivotCol) {
933  assert(pivotCol >= getNumFixedCols() && "Refusing to pivot invalid column");
934  assert(!unknownFromColumn(pivotCol).isSymbol);
935 
936  swapRowWithCol(pivotRow, pivotCol);
937  std::swap(tableau(pivotRow, 0), tableau(pivotRow, pivotCol));
938  // We need to negate the whole pivot row except for the pivot column.
939  if (tableau(pivotRow, 0) < 0) {
940  // If the denominator is negative, we negate the row by simply negating the
941  // denominator.
942  tableau(pivotRow, 0) = -tableau(pivotRow, 0);
943  tableau(pivotRow, pivotCol) = -tableau(pivotRow, pivotCol);
944  } else {
945  for (unsigned col = 1, e = getNumColumns(); col < e; ++col) {
946  if (col == pivotCol)
947  continue;
948  tableau(pivotRow, col) = -tableau(pivotRow, col);
949  }
950  }
951  tableau.normalizeRow(pivotRow);
952 
953  for (unsigned row = 0, numRows = getNumRows(); row < numRows; ++row) {
954  if (row == pivotRow)
955  continue;
956  if (tableau(row, pivotCol) == 0) // Nothing to do.
957  continue;
958  tableau(row, 0) *= tableau(pivotRow, 0);
959  for (unsigned col = 1, numCols = getNumColumns(); col < numCols; ++col) {
960  if (col == pivotCol)
961  continue;
962  // Add rather than subtract because the pivot row has been negated.
963  tableau(row, col) = tableau(row, col) * tableau(pivotRow, 0) +
964  tableau(row, pivotCol) * tableau(pivotRow, col);
965  }
966  tableau(row, pivotCol) *= tableau(pivotRow, pivotCol);
967  tableau.normalizeRow(row);
968  }
969 }
970 
971 /// Perform pivots until the unknown has a non-negative sample value or until
972 /// no more upward pivots can be performed. Return success if we were able to
973 /// bring the row to a non-negative sample value, and failure otherwise.
974 LogicalResult Simplex::restoreRow(Unknown &u) {
975  assert(u.orientation == Orientation::Row &&
976  "unknown should be in row position");
977 
978  while (tableau(u.pos, 1) < 0) {
979  Optional<Pivot> maybePivot = findPivot(u.pos, Direction::Up);
980  if (!maybePivot)
981  break;
982 
983  pivot(*maybePivot);
984  if (u.orientation == Orientation::Column)
985  return success(); // the unknown is unbounded above.
986  }
987  return success(tableau(u.pos, 1) >= 0);
988 }
989 
990 /// Find a row that can be used to pivot the column in the specified direction.
991 /// This returns an empty optional if and only if the column is unbounded in the
992 /// specified direction (ignoring skipRow, if skipRow is set).
993 ///
994 /// If skipRow is set, this row is not considered, and (if it is restricted) its
995 /// restriction may be violated by the returned pivot. Usually, skipRow is set
996 /// because we don't want to move it to column position unless it is unbounded,
997 /// and we are either trying to increase the value of skipRow or explicitly
998 /// trying to make skipRow negative, so we are not concerned about this.
999 ///
1000 /// If the direction is up (resp. down) and a restricted row has a negative
1001 /// (positive) coefficient for the column, then this row imposes a bound on how
1002 /// much the sample value of the column can change. Such a row with constant
1003 /// term c and coefficient f for the column imposes a bound of c/|f| on the
1004 /// change in sample value (in the specified direction). (note that c is
1005 /// non-negative here since the row is restricted and the tableau is consistent)
1006 ///
1007 /// We iterate through the rows and pick the row which imposes the most
1008 /// stringent bound, since pivoting with a row changes the row's sample value to
1009 /// 0 and hence saturates the bound it imposes. We break ties between rows that
1010 /// impose the same bound by considering a lexicographic ordering where we
1011 /// prefer unknowns with lower index value.
1012 Optional<unsigned> Simplex::findPivotRow(Optional<unsigned> skipRow,
1013  Direction direction,
1014  unsigned col) const {
1015  Optional<unsigned> retRow;
1016  // Initialize these to zero in order to silence a warning about retElem and
1017  // retConst being used uninitialized in the initialization of `diff` below. In
1018  // reality, these are always initialized when that line is reached since these
1019  // are set whenever retRow is set.
1020  MPInt retElem, retConst;
1021  for (unsigned row = nRedundant, e = getNumRows(); row < e; ++row) {
1022  if (skipRow && row == *skipRow)
1023  continue;
1024  MPInt elem = tableau(row, col);
1025  if (elem == 0)
1026  continue;
1027  if (!unknownFromRow(row).restricted)
1028  continue;
1029  if (signMatchesDirection(elem, direction))
1030  continue;
1031  MPInt constTerm = tableau(row, 1);
1032 
1033  if (!retRow) {
1034  retRow = row;
1035  retElem = elem;
1036  retConst = constTerm;
1037  continue;
1038  }
1039 
1040  MPInt diff = retConst * elem - constTerm * retElem;
1041  if ((diff == 0 && rowUnknown[row] < rowUnknown[*retRow]) ||
1042  (diff != 0 && !signMatchesDirection(diff, direction))) {
1043  retRow = row;
1044  retElem = elem;
1045  retConst = constTerm;
1046  }
1047  }
1048  return retRow;
1049 }
1050 
1051 bool SimplexBase::isEmpty() const { return empty; }
1052 
1053 void SimplexBase::swapRows(unsigned i, unsigned j) {
1054  if (i == j)
1055  return;
1056  tableau.swapRows(i, j);
1057  std::swap(rowUnknown[i], rowUnknown[j]);
1058  unknownFromRow(i).pos = i;
1059  unknownFromRow(j).pos = j;
1060 }
1061 
1062 void SimplexBase::swapColumns(unsigned i, unsigned j) {
1063  assert(i < getNumColumns() && j < getNumColumns() &&
1064  "Invalid columns provided!");
1065  if (i == j)
1066  return;
1067  tableau.swapColumns(i, j);
1068  std::swap(colUnknown[i], colUnknown[j]);
1069  unknownFromColumn(i).pos = i;
1070  unknownFromColumn(j).pos = j;
1071 }
1072 
1073 /// Mark this tableau empty and push an entry to the undo stack.
1075  // If the set is already empty, then we shouldn't add another UnmarkEmpty log
1076  // entry, since in that case the Simplex will be erroneously marked as
1077  // non-empty when rolling back past this point.
1078  if (empty)
1079  return;
1081  empty = true;
1082 }
1083 
1084 /// Add an inequality to the tableau. If coeffs is c_0, c_1, ... c_n, where n
1085 /// is the current number of variables, then the corresponding inequality is
1086 /// c_n + c_0*x_0 + c_1*x_1 + ... + c_{n-1}*x_{n-1} >= 0.
1087 ///
1088 /// We add the inequality and mark it as restricted. We then try to make its
1089 /// sample value non-negative. If this is not possible, the tableau has become
1090 /// empty and we mark it as such.
1092  unsigned conIndex = addRow(coeffs, /*makeRestricted=*/true);
1093  LogicalResult result = restoreRow(con[conIndex]);
1094  if (failed(result))
1095  markEmpty();
1096 }
1097 
1098 /// Add an equality to the tableau. If coeffs is c_0, c_1, ... c_n, where n
1099 /// is the current number of variables, then the corresponding equality is
1100 /// c_n + c_0*x_0 + c_1*x_1 + ... + c_{n-1}*x_{n-1} == 0.
1101 ///
1102 /// We simply add two opposing inequalities, which force the expression to
1103 /// be zero.
1105  addInequality(coeffs);
1106  SmallVector<MPInt, 8> negatedCoeffs;
1107  for (const MPInt &coeff : coeffs)
1108  negatedCoeffs.emplace_back(-coeff);
1109  addInequality(negatedCoeffs);
1110 }
1111 
1112 unsigned SimplexBase::getNumVariables() const { return var.size(); }
1113 unsigned SimplexBase::getNumConstraints() const { return con.size(); }
1114 
1115 /// Return a snapshot of the current state. This is just the current size of the
1116 /// undo log.
1117 unsigned SimplexBase::getSnapshot() const { return undoLog.size(); }
1118 
1120  SmallVector<int, 8> basis;
1121  for (int index : colUnknown) {
1122  if (index != nullIndex)
1123  basis.push_back(index);
1124  }
1125  savedBases.push_back(std::move(basis));
1126 
1127  undoLog.emplace_back(UndoLogEntry::RestoreBasis);
1128  return undoLog.size() - 1;
1129 }
1130 
1132  assert(con.back().orientation == Orientation::Row);
1133 
1134  // Move this unknown to the last row and remove the last row from the
1135  // tableau.
1136  swapRows(con.back().pos, getNumRows() - 1);
1137  // It is not strictly necessary to shrink the tableau, but for now we
1138  // maintain the invariant that the tableau has exactly getNumRows()
1139  // rows.
1141  rowUnknown.pop_back();
1142  con.pop_back();
1143 }
1144 
1145 // This doesn't find a pivot row only if the column has zero
1146 // coefficients for every row.
1147 //
1148 // If the unknown is a constraint, this can't happen, since it was added
1149 // initially as a row. Such a row could never have been pivoted to a column. So
1150 // a pivot row will always be found if we have a constraint.
1151 //
1152 // If we have a variable, then the column has zero coefficients for every row
1153 // iff no constraints have been added with a non-zero coefficient for this row.
1155  for (unsigned row = nRedundant, e = getNumRows(); row < e; ++row)
1156  if (tableau(row, col) != 0)
1157  return row;
1158  return {};
1159 }
1160 
1161 // It's not valid to remove the constraint by deleting the column since this
1162 // would result in an invalid basis.
1163 void Simplex::undoLastConstraint() {
1164  if (con.back().orientation == Orientation::Column) {
1165  // We try to find any pivot row for this column that preserves tableau
1166  // consistency (except possibly the column itself, which is going to be
1167  // deallocated anyway).
1168  //
1169  // If no pivot row is found in either direction, then the unknown is
1170  // unbounded in both directions and we are free to perform any pivot at
1171  // all. To do this, we just need to find any row with a non-zero
1172  // coefficient for the column. findAnyPivotRow will always be able to
1173  // find such a row for a constraint.
1174  unsigned column = con.back().pos;
1175  if (Optional<unsigned> maybeRow = findPivotRow({}, Direction::Up, column)) {
1176  pivot(*maybeRow, column);
1177  } else if (Optional<unsigned> maybeRow =
1178  findPivotRow({}, Direction::Down, column)) {
1179  pivot(*maybeRow, column);
1180  } else {
1181  Optional<unsigned> row = findAnyPivotRow(column);
1182  assert(row && "Pivot should always exist for a constraint!");
1183  pivot(*row, column);
1184  }
1185  }
1187 }
1188 
1189 // It's not valid to remove the constraint by deleting the column since this
1190 // would result in an invalid basis.
1192  if (con.back().orientation == Orientation::Column) {
1193  // When removing the last constraint during a rollback, we just need to find
1194  // any pivot at all, i.e., any row with non-zero coefficient for the
1195  // column, because when rolling back a lexicographic simplex, we always
1196  // end by restoring the exact basis that was present at the time of the
1197  // snapshot, so what pivots we perform while undoing doesn't matter as
1198  // long as we get the unknown to row orientation and remove it.
1199  unsigned column = con.back().pos;
1200  Optional<unsigned> row = findAnyPivotRow(column);
1201  assert(row && "Pivot should always exist for a constraint!");
1202  pivot(*row, column);
1203  }
1205 }
1206 
1208  if (entry == UndoLogEntry::RemoveLastConstraint) {
1209  // Simplex and LexSimplex handle this differently, so we call out to a
1210  // virtual function to handle this.
1212  } else if (entry == UndoLogEntry::RemoveLastVariable) {
1213  // Whenever we are rolling back the addition of a variable, it is guaranteed
1214  // that the variable will be in column position.
1215  //
1216  // We can see this as follows: any constraint that depends on this variable
1217  // was added after this variable was added, so the addition of such
1218  // constraints should already have been rolled back by the time we get to
1219  // rolling back the addition of the variable. Therefore, no constraint
1220  // currently has a component along the variable, so the variable itself must
1221  // be part of the basis.
1222  assert(var.back().orientation == Orientation::Column &&
1223  "Variable to be removed must be in column orientation!");
1224 
1225  if (var.back().isSymbol)
1226  nSymbol--;
1227 
1228  // Move this variable to the last column and remove the column from the
1229  // tableau.
1230  swapColumns(var.back().pos, getNumColumns() - 1);
1232  var.pop_back();
1233  colUnknown.pop_back();
1234  } else if (entry == UndoLogEntry::UnmarkEmpty) {
1235  empty = false;
1236  } else if (entry == UndoLogEntry::UnmarkLastRedundant) {
1237  nRedundant--;
1238  } else if (entry == UndoLogEntry::RestoreBasis) {
1239  assert(!savedBases.empty() && "No bases saved!");
1240 
1241  SmallVector<int, 8> basis = std::move(savedBases.back());
1242  savedBases.pop_back();
1243 
1244  for (int index : basis) {
1245  Unknown &u = unknownFromIndex(index);
1247  continue;
1248  for (unsigned col = getNumFixedCols(), e = getNumColumns(); col < e;
1249  col++) {
1250  assert(colUnknown[col] != nullIndex &&
1251  "Column should not be a fixed column!");
1252  if (llvm::is_contained(basis, colUnknown[col]))
1253  continue;
1254  if (tableau(u.pos, col) == 0)
1255  continue;
1256  pivot(u.pos, col);
1257  break;
1258  }
1259 
1260  assert(u.orientation == Orientation::Column && "No pivot found!");
1261  }
1262  }
1263 }
1264 
1265 /// Rollback to the specified snapshot.
1266 ///
1267 /// We undo all the log entries until the log size when the snapshot was taken
1268 /// is reached.
1269 void SimplexBase::rollback(unsigned snapshot) {
1270  while (undoLog.size() > snapshot) {
1271  undo(undoLog.back());
1272  undoLog.pop_back();
1273  }
1274 }
1275 
1276 /// We add the usual floor division constraints:
1277 /// `0 <= coeffs - denom*q <= denom - 1`, where `q` is the new division
1278 /// variable.
1279 ///
1280 /// This constrains the remainder `coeffs - denom*q` to be in the
1281 /// range `[0, denom - 1]`, which fixes the integer value of the quotient `q`.
1283  const MPInt &denom) {
1284  assert(denom > 0 && "Denominator must be positive!");
1285  appendVariable();
1286 
1287  SmallVector<MPInt, 8> ineq(coeffs.begin(), coeffs.end());
1288  MPInt constTerm = ineq.back();
1289  ineq.back() = -denom;
1290  ineq.push_back(constTerm);
1291  addInequality(ineq);
1292 
1293  for (MPInt &coeff : ineq)
1294  coeff = -coeff;
1295  ineq.back() += denom - 1;
1296  addInequality(ineq);
1297 }
1298 
1299 void SimplexBase::appendVariable(unsigned count) {
1300  if (count == 0)
1301  return;
1302  var.reserve(var.size() + count);
1303  colUnknown.reserve(colUnknown.size() + count);
1304  for (unsigned i = 0; i < count; ++i) {
1305  var.emplace_back(Orientation::Column, /*restricted=*/false,
1306  /*pos=*/getNumColumns() + i);
1307  colUnknown.push_back(var.size() - 1);
1308  }
1310  undoLog.insert(undoLog.end(), count, UndoLogEntry::RemoveLastVariable);
1311 }
1312 
1313 /// Add all the constraints from the given IntegerRelation.
1315  assert(rel.getNumVars() == getNumVariables() &&
1316  "IntegerRelation must have same dimensionality as simplex");
1317  for (unsigned i = 0, e = rel.getNumInequalities(); i < e; ++i)
1318  addInequality(rel.getInequality(i));
1319  for (unsigned i = 0, e = rel.getNumEqualities(); i < e; ++i)
1320  addEquality(rel.getEquality(i));
1321 }
1322 
1324  unsigned row) {
1325  // Keep trying to find a pivot for the row in the specified direction.
1326  while (Optional<Pivot> maybePivot = findPivot(row, direction)) {
1327  // If findPivot returns a pivot involving the row itself, then the optimum
1328  // is unbounded, so we return std::nullopt.
1329  if (maybePivot->row == row)
1330  return OptimumKind::Unbounded;
1331  pivot(*maybePivot);
1332  }
1333 
1334  // The row has reached its optimal sample value, which we return.
1335  // The sample value is the entry in the constant column divided by the common
1336  // denominator for this row.
1337  return Fraction(tableau(row, 1), tableau(row, 0));
1338 }
1339 
1340 /// Compute the optimum of the specified expression in the specified direction,
1341 /// or std::nullopt if it is unbounded.
1343  ArrayRef<MPInt> coeffs) {
1344  if (empty)
1345  return OptimumKind::Empty;
1346 
1347  SimplexRollbackScopeExit scopeExit(*this);
1348  unsigned conIndex = addRow(coeffs);
1349  unsigned row = con[conIndex].pos;
1350  return computeRowOptimum(direction, row);
1351 }
1352 
1354  Unknown &u) {
1355  if (empty)
1356  return OptimumKind::Empty;
1357  if (u.orientation == Orientation::Column) {
1358  unsigned column = u.pos;
1359  Optional<unsigned> pivotRow = findPivotRow({}, direction, column);
1360  // If no pivot is returned, the constraint is unbounded in the specified
1361  // direction.
1362  if (!pivotRow)
1363  return OptimumKind::Unbounded;
1364  pivot(*pivotRow, column);
1365  }
1366 
1367  unsigned row = u.pos;
1368  MaybeOptimum<Fraction> optimum = computeRowOptimum(direction, row);
1369  if (u.restricted && direction == Direction::Down &&
1370  (optimum.isUnbounded() || *optimum < Fraction(0, 1))) {
1371  if (failed(restoreRow(u)))
1372  llvm_unreachable("Could not restore row!");
1373  }
1374  return optimum;
1375 }
1376 
1377 bool Simplex::isBoundedAlongConstraint(unsigned constraintIndex) {
1378  assert(!empty && "It is not meaningful to ask whether a direction is bounded "
1379  "in an empty set.");
1380  // The constraint's perpendicular is already bounded below, since it is a
1381  // constraint. If it is also bounded above, we can return true.
1382  return computeOptimum(Direction::Up, con[constraintIndex]).isBounded();
1383 }
1384 
1385 /// Redundant constraints are those that are in row orientation and lie in
1386 /// rows 0 to nRedundant - 1.
1387 bool Simplex::isMarkedRedundant(unsigned constraintIndex) const {
1388  const Unknown &u = con[constraintIndex];
1389  return u.orientation == Orientation::Row && u.pos < nRedundant;
1390 }
1391 
1392 /// Mark the specified row redundant.
1393 ///
1394 /// This is done by moving the unknown to the end of the block of redundant
1395 /// rows (namely, to row nRedundant) and incrementing nRedundant to
1396 /// accomodate the new redundant row.
1397 void Simplex::markRowRedundant(Unknown &u) {
1398  assert(u.orientation == Orientation::Row &&
1399  "Unknown should be in row position!");
1400  assert(u.pos >= nRedundant && "Unknown is already marked redundant!");
1401  swapRows(u.pos, nRedundant);
1402  ++nRedundant;
1404 }
1405 
1406 /// Find a subset of constraints that is redundant and mark them redundant.
1407 void Simplex::detectRedundant(unsigned offset, unsigned count) {
1408  assert(offset + count <= con.size() && "invalid range!");
1409  // It is not meaningful to talk about redundancy for empty sets.
1410  if (empty)
1411  return;
1412 
1413  // Iterate through the constraints and check for each one if it can attain
1414  // negative sample values. If it can, it's not redundant. Otherwise, it is.
1415  // We mark redundant constraints redundant.
1416  //
1417  // Constraints that get marked redundant in one iteration are not respected
1418  // when checking constraints in later iterations. This prevents, for example,
1419  // two identical constraints both being marked redundant since each is
1420  // redundant given the other one. In this example, only the first of the
1421  // constraints that is processed will get marked redundant, as it should be.
1422  for (unsigned i = 0; i < count; ++i) {
1423  Unknown &u = con[offset + i];
1424  if (u.orientation == Orientation::Column) {
1425  unsigned column = u.pos;
1426  Optional<unsigned> pivotRow = findPivotRow({}, Direction::Down, column);
1427  // If no downward pivot is returned, the constraint is unbounded below
1428  // and hence not redundant.
1429  if (!pivotRow)
1430  continue;
1431  pivot(*pivotRow, column);
1432  }
1433 
1434  unsigned row = u.pos;
1436  if (minimum.isUnbounded() || *minimum < Fraction(0, 1)) {
1437  // Constraint is unbounded below or can attain negative sample values and
1438  // hence is not redundant.
1439  if (failed(restoreRow(u)))
1440  llvm_unreachable("Could not restore non-redundant row!");
1441  continue;
1442  }
1443 
1444  markRowRedundant(u);
1445  }
1446 }
1447 
1449  if (empty)
1450  return false;
1451 
1452  SmallVector<MPInt, 8> dir(var.size() + 1);
1453  for (unsigned i = 0; i < var.size(); ++i) {
1454  dir[i] = 1;
1455 
1457  return true;
1458 
1460  return true;
1461 
1462  dir[i] = 0;
1463  }
1464  return false;
1465 }
1466 
1467 /// Make a tableau to represent a pair of points in the original tableau.
1468 ///
1469 /// The product constraints and variables are stored as: first A's, then B's.
1470 ///
1471 /// The product tableau has row layout:
1472 /// A's redundant rows, B's redundant rows, A's other rows, B's other rows.
1473 ///
1474 /// It has column layout:
1475 /// denominator, constant, A's columns, B's columns.
1477  unsigned numVar = a.getNumVariables() + b.getNumVariables();
1478  unsigned numCon = a.getNumConstraints() + b.getNumConstraints();
1479  Simplex result(numVar);
1480 
1481  result.tableau.reserveRows(numCon);
1482  result.empty = a.empty || b.empty;
1483 
1484  auto concat = [](ArrayRef<Unknown> v, ArrayRef<Unknown> w) {
1485  SmallVector<Unknown, 8> result;
1486  result.reserve(v.size() + w.size());
1487  result.insert(result.end(), v.begin(), v.end());
1488  result.insert(result.end(), w.begin(), w.end());
1489  return result;
1490  };
1491  result.con = concat(a.con, b.con);
1492  result.var = concat(a.var, b.var);
1493 
1494  auto indexFromBIndex = [&](int index) {
1495  return index >= 0 ? a.getNumVariables() + index
1496  : ~(a.getNumConstraints() + ~index);
1497  };
1498 
1499  result.colUnknown.assign(2, nullIndex);
1500  for (unsigned i = 2, e = a.getNumColumns(); i < e; ++i) {
1501  result.colUnknown.push_back(a.colUnknown[i]);
1502  result.unknownFromIndex(result.colUnknown.back()).pos =
1503  result.colUnknown.size() - 1;
1504  }
1505  for (unsigned i = 2, e = b.getNumColumns(); i < e; ++i) {
1506  result.colUnknown.push_back(indexFromBIndex(b.colUnknown[i]));
1507  result.unknownFromIndex(result.colUnknown.back()).pos =
1508  result.colUnknown.size() - 1;
1509  }
1510 
1511  auto appendRowFromA = [&](unsigned row) {
1512  unsigned resultRow = result.tableau.appendExtraRow();
1513  for (unsigned col = 0, e = a.getNumColumns(); col < e; ++col)
1514  result.tableau(resultRow, col) = a.tableau(row, col);
1515  result.rowUnknown.push_back(a.rowUnknown[row]);
1516  result.unknownFromIndex(result.rowUnknown.back()).pos =
1517  result.rowUnknown.size() - 1;
1518  };
1519 
1520  // Also fixes the corresponding entry in rowUnknown and var/con (as the case
1521  // may be).
1522  auto appendRowFromB = [&](unsigned row) {
1523  unsigned resultRow = result.tableau.appendExtraRow();
1524  result.tableau(resultRow, 0) = b.tableau(row, 0);
1525  result.tableau(resultRow, 1) = b.tableau(row, 1);
1526 
1527  unsigned offset = a.getNumColumns() - 2;
1528  for (unsigned col = 2, e = b.getNumColumns(); col < e; ++col)
1529  result.tableau(resultRow, offset + col) = b.tableau(row, col);
1530  result.rowUnknown.push_back(indexFromBIndex(b.rowUnknown[row]));
1531  result.unknownFromIndex(result.rowUnknown.back()).pos =
1532  result.rowUnknown.size() - 1;
1533  };
1534 
1535  result.nRedundant = a.nRedundant + b.nRedundant;
1536  for (unsigned row = 0; row < a.nRedundant; ++row)
1537  appendRowFromA(row);
1538  for (unsigned row = 0; row < b.nRedundant; ++row)
1539  appendRowFromB(row);
1540  for (unsigned row = a.nRedundant, e = a.getNumRows(); row < e; ++row)
1541  appendRowFromA(row);
1542  for (unsigned row = b.nRedundant, e = b.getNumRows(); row < e; ++row)
1543  appendRowFromB(row);
1544 
1545  return result;
1546 }
1547 
1549  if (empty)
1550  return {};
1551 
1552  SmallVector<Fraction, 8> sample;
1553  sample.reserve(var.size());
1554  // Push the sample value for each variable into the vector.
1555  for (const Unknown &u : var) {
1556  if (u.orientation == Orientation::Column) {
1557  // If the variable is in column position, its sample value is zero.
1558  sample.emplace_back(0, 1);
1559  } else {
1560  // If the variable is in row position, its sample value is the
1561  // entry in the constant column divided by the denominator.
1562  MPInt denom = tableau(u.pos, 0);
1563  sample.emplace_back(tableau(u.pos, 1), denom);
1564  }
1565  }
1566  return sample;
1567 }
1568 
1570  addRow(coeffs, /*makeRestricted=*/true);
1571 }
1572 
1573 MaybeOptimum<SmallVector<Fraction, 8>> LexSimplex::getRationalSample() const {
1574  if (empty)
1575  return OptimumKind::Empty;
1576 
1577  SmallVector<Fraction, 8> sample;
1578  sample.reserve(var.size());
1579  // Push the sample value for each variable into the vector.
1580  for (const Unknown &u : var) {
1581  // When the big M parameter is being used, each variable x is represented
1582  // as M + x, so its sample value is finite if and only if it is of the
1583  // form 1*M + c. If the coefficient of M is not one then the sample value
1584  // is infinite, and we return an empty optional.
1585 
1586  if (u.orientation == Orientation::Column) {
1587  // If the variable is in column position, the sample value of M + x is
1588  // zero, so x = -M which is unbounded.
1589  return OptimumKind::Unbounded;
1590  }
1591 
1592  // If the variable is in row position, its sample value is the
1593  // entry in the constant column divided by the denominator.
1594  MPInt denom = tableau(u.pos, 0);
1595  if (usingBigM)
1596  if (tableau(u.pos, 2) != denom)
1597  return OptimumKind::Unbounded;
1598  sample.emplace_back(tableau(u.pos, 1), denom);
1599  }
1600  return sample;
1601 }
1602 
1604  // If the tableau is empty, no sample point exists.
1605  if (empty)
1606  return {};
1607 
1608  // The value will always exist since the Simplex is non-empty.
1609  SmallVector<Fraction, 8> rationalSample = *getRationalSample();
1610  SmallVector<MPInt, 8> integerSample;
1611  integerSample.reserve(var.size());
1612  for (const Fraction &coord : rationalSample) {
1613  // If the sample is non-integral, return std::nullopt.
1614  if (coord.num % coord.den != 0)
1615  return {};
1616  integerSample.push_back(coord.num / coord.den);
1617  }
1618  return integerSample;
1619 }
1620 
1621 /// Given a simplex for a polytope, construct a new simplex whose variables are
1622 /// identified with a pair of points (x, y) in the original polytope. Supports
1623 /// some operations needed for generalized basis reduction. In what follows,
1624 /// dotProduct(x, y) = x_1 * y_1 + x_2 * y_2 + ... x_n * y_n where n is the
1625 /// dimension of the original polytope.
1626 ///
1627 /// This supports adding equality constraints dotProduct(dir, x - y) == 0. It
1628 /// also supports rolling back this addition, by maintaining a snapshot stack
1629 /// that contains a snapshot of the Simplex's state for each equality, just
1630 /// before that equality was added.
1633 
1634 public:
1635  GBRSimplex(const Simplex &originalSimplex)
1636  : simplex(Simplex::makeProduct(originalSimplex, originalSimplex)),
1637  simplexConstraintOffset(simplex.getNumConstraints()) {}
1638 
1639  /// Add an equality dotProduct(dir, x - y) == 0.
1640  /// First pushes a snapshot for the current simplex state to the stack so
1641  /// that this can be rolled back later.
1643  assert(llvm::any_of(dir, [](const MPInt &x) { return x != 0; }) &&
1644  "Direction passed is the zero vector!");
1645  snapshotStack.push_back(simplex.getSnapshot());
1646  simplex.addEquality(getCoeffsForDirection(dir));
1647  }
1648  /// Compute max(dotProduct(dir, x - y)).
1650  MaybeOptimum<Fraction> maybeWidth =
1651  simplex.computeOptimum(Direction::Up, getCoeffsForDirection(dir));
1652  assert(maybeWidth.isBounded() && "Width should be bounded!");
1653  return *maybeWidth;
1654  }
1655 
1656  /// Compute max(dotProduct(dir, x - y)) and save the dual variables for only
1657  /// the direction equalities to `dual`.
1659  SmallVectorImpl<MPInt> &dual,
1660  MPInt &dualDenom) {
1661  // We can't just call into computeWidth or computeOptimum since we need to
1662  // access the state of the tableau after computing the optimum, and these
1663  // functions rollback the insertion of the objective function into the
1664  // tableau before returning. We instead add a row for the objective function
1665  // ourselves, call into computeOptimum, compute the duals from the tableau
1666  // state, and finally rollback the addition of the row before returning.
1667  SimplexRollbackScopeExit scopeExit(simplex);
1668  unsigned conIndex = simplex.addRow(getCoeffsForDirection(dir));
1669  unsigned row = simplex.con[conIndex].pos;
1670  MaybeOptimum<Fraction> maybeWidth =
1671  simplex.computeRowOptimum(Simplex::Direction::Up, row);
1672  assert(maybeWidth.isBounded() && "Width should be bounded!");
1673  dualDenom = simplex.tableau(row, 0);
1674  dual.clear();
1675 
1676  // The increment is i += 2 because equalities are added as two inequalities,
1677  // one positive and one negative. Each iteration processes one equality.
1678  for (unsigned i = simplexConstraintOffset; i < conIndex; i += 2) {
1679  // The dual variable for an inequality in column orientation is the
1680  // negative of its coefficient at the objective row. If the inequality is
1681  // in row orientation, the corresponding dual variable is zero.
1682  //
1683  // We want the dual for the original equality, which corresponds to two
1684  // inequalities: a positive inequality, which has the same coefficients as
1685  // the equality, and a negative equality, which has negated coefficients.
1686  //
1687  // Note that at most one of these inequalities can be in column
1688  // orientation because the column unknowns should form a basis and hence
1689  // must be linearly independent. If the positive inequality is in column
1690  // position, its dual is the dual corresponding to the equality. If the
1691  // negative inequality is in column position, the negation of its dual is
1692  // the dual corresponding to the equality. If neither is in column
1693  // position, then that means that this equality is redundant, and its dual
1694  // is zero.
1695  //
1696  // Note that it is NOT valid to perform pivots during the computation of
1697  // the duals. This entire dual computation must be performed on the same
1698  // tableau configuration.
1699  assert(!(simplex.con[i].orientation == Orientation::Column &&
1700  simplex.con[i + 1].orientation == Orientation::Column) &&
1701  "Both inequalities for the equality cannot be in column "
1702  "orientation!");
1703  if (simplex.con[i].orientation == Orientation::Column)
1704  dual.push_back(-simplex.tableau(row, simplex.con[i].pos));
1705  else if (simplex.con[i + 1].orientation == Orientation::Column)
1706  dual.push_back(simplex.tableau(row, simplex.con[i + 1].pos));
1707  else
1708  dual.emplace_back(0);
1709  }
1710  return *maybeWidth;
1711  }
1712 
1713  /// Remove the last equality that was added through addEqualityForDirection.
1714  ///
1715  /// We do this by rolling back to the snapshot at the top of the stack, which
1716  /// should be a snapshot taken just before the last equality was added.
1718  assert(!snapshotStack.empty() && "Snapshot stack is empty!");
1719  simplex.rollback(snapshotStack.back());
1720  snapshotStack.pop_back();
1721  }
1722 
1723 private:
1724  /// Returns coefficients of the expression 'dot_product(dir, x - y)',
1725  /// i.e., dir_1 * x_1 + dir_2 * x_2 + ... + dir_n * x_n
1726  /// - dir_1 * y_1 - dir_2 * y_2 - ... - dir_n * y_n,
1727  /// where n is the dimension of the original polytope.
1728  SmallVector<MPInt, 8> getCoeffsForDirection(ArrayRef<MPInt> dir) {
1729  assert(2 * dir.size() == simplex.getNumVariables() &&
1730  "Direction vector has wrong dimensionality");
1731  SmallVector<MPInt, 8> coeffs(dir.begin(), dir.end());
1732  coeffs.reserve(2 * dir.size());
1733  for (const MPInt &coeff : dir)
1734  coeffs.push_back(-coeff);
1735  coeffs.emplace_back(0); // constant term
1736  return coeffs;
1737  }
1738 
1739  Simplex simplex;
1740  /// The first index of the equality constraints, the index immediately after
1741  /// the last constraint in the initial product simplex.
1742  unsigned simplexConstraintOffset;
1743  /// A stack of snapshots, used for rolling back.
1744  SmallVector<unsigned, 8> snapshotStack;
1745 };
1746 
1747 /// Reduce the basis to try and find a direction in which the polytope is
1748 /// "thin". This only works for bounded polytopes.
1749 ///
1750 /// This is an implementation of the algorithm described in the paper
1751 /// "An Implementation of Generalized Basis Reduction for Integer Programming"
1752 /// by W. Cook, T. Rutherford, H. E. Scarf, D. Shallcross.
1753 ///
1754 /// Let b_{level}, b_{level + 1}, ... b_n be the current basis.
1755 /// Let width_i(v) = max <v, x - y> where x and y are points in the original
1756 /// polytope such that <b_j, x - y> = 0 is satisfied for all level <= j < i.
1757 ///
1758 /// In every iteration, we first replace b_{i+1} with b_{i+1} + u*b_i, where u
1759 /// is the integer such that width_i(b_{i+1} + u*b_i) is minimized. Let dual_i
1760 /// be the dual variable associated with the constraint <b_i, x - y> = 0 when
1761 /// computing width_{i+1}(b_{i+1}). It can be shown that dual_i is the
1762 /// minimizing value of u, if it were allowed to be fractional. Due to
1763 /// convexity, the minimizing integer value is either floor(dual_i) or
1764 /// ceil(dual_i), so we just need to check which of these gives a lower
1765 /// width_{i+1} value. If dual_i turned out to be an integer, then u = dual_i.
1766 ///
1767 /// Now if width_i(b_{i+1}) < 0.75 * width_i(b_i), we swap b_i and (the new)
1768 /// b_{i + 1} and decrement i (unless i = level, in which case we stay at the
1769 /// same i). Otherwise, we increment i.
1770 ///
1771 /// We keep f values and duals cached and invalidate them when necessary.
1772 /// Whenever possible, we use them instead of recomputing them. We implement the
1773 /// algorithm as follows.
1774 ///
1775 /// In an iteration at i we need to compute:
1776 /// a) width_i(b_{i + 1})
1777 /// b) width_i(b_i)
1778 /// c) the integer u that minimizes width_i(b_{i + 1} + u*b_i)
1779 ///
1780 /// If width_i(b_i) is not already cached, we compute it.
1781 ///
1782 /// If the duals are not already cached, we compute width_{i+1}(b_{i+1}) and
1783 /// store the duals from this computation.
1784 ///
1785 /// We call updateBasisWithUAndGetFCandidate, which finds the minimizing value
1786 /// of u as explained before, caches the duals from this computation, sets
1787 /// b_{i+1} to b_{i+1} + u*b_i, and returns the new value of width_i(b_{i+1}).
1788 ///
1789 /// Now if width_i(b_{i+1}) < 0.75 * width_i(b_i), we swap b_i and b_{i+1} and
1790 /// decrement i, resulting in the basis
1791 /// ... b_{i - 1}, b_{i + 1} + u*b_i, b_i, b_{i+2}, ...
1792 /// with corresponding f values
1793 /// ... width_{i-1}(b_{i-1}), width_i(b_{i+1} + u*b_i), width_{i+1}(b_i), ...
1794 /// The values up to i - 1 remain unchanged. We have just gotten the middle
1795 /// value from updateBasisWithUAndGetFCandidate, so we can update that in the
1796 /// cache. The value at width_{i+1}(b_i) is unknown, so we evict this value from
1797 /// the cache. The iteration after decrementing needs exactly the duals from the
1798 /// computation of width_i(b_{i + 1} + u*b_i), so we keep these in the cache.
1799 ///
1800 /// When incrementing i, no cached f values get invalidated. However, the cached
1801 /// duals do get invalidated as the duals for the higher levels are different.
1802 void Simplex::reduceBasis(Matrix &basis, unsigned level) {
1803  const Fraction epsilon(3, 4);
1804 
1805  if (level == basis.getNumRows() - 1)
1806  return;
1807 
1808  GBRSimplex gbrSimplex(*this);
1810  SmallVector<MPInt, 8> dual;
1811  MPInt dualDenom;
1812 
1813  // Finds the value of u that minimizes width_i(b_{i+1} + u*b_i), caches the
1814  // duals from this computation, sets b_{i+1} to b_{i+1} + u*b_i, and returns
1815  // the new value of width_i(b_{i+1}).
1816  //
1817  // If dual_i is not an integer, the minimizing value must be either
1818  // floor(dual_i) or ceil(dual_i). We compute the expression for both and
1819  // choose the minimizing value.
1820  //
1821  // If dual_i is an integer, we don't need to perform these computations. We
1822  // know that in this case,
1823  // a) u = dual_i.
1824  // b) one can show that dual_j for j < i are the same duals we would have
1825  // gotten from computing width_i(b_{i + 1} + u*b_i), so the correct duals
1826  // are the ones already in the cache.
1827  // c) width_i(b_{i+1} + u*b_i) = min_{alpha} width_i(b_{i+1} + alpha * b_i),
1828  // which
1829  // one can show is equal to width_{i+1}(b_{i+1}). The latter value must
1830  // be in the cache, so we get it from there and return it.
1831  auto updateBasisWithUAndGetFCandidate = [&](unsigned i) -> Fraction {
1832  assert(i < level + dual.size() && "dual_i is not known!");
1833 
1834  MPInt u = floorDiv(dual[i - level], dualDenom);
1835  basis.addToRow(i, i + 1, u);
1836  if (dual[i - level] % dualDenom != 0) {
1837  SmallVector<MPInt, 8> candidateDual[2];
1838  MPInt candidateDualDenom[2];
1839  Fraction widthI[2];
1840 
1841  // Initially u is floor(dual) and basis reflects this.
1842  widthI[0] = gbrSimplex.computeWidthAndDuals(
1843  basis.getRow(i + 1), candidateDual[0], candidateDualDenom[0]);
1844 
1845  // Now try ceil(dual), i.e. floor(dual) + 1.
1846  ++u;
1847  basis.addToRow(i, i + 1, 1);
1848  widthI[1] = gbrSimplex.computeWidthAndDuals(
1849  basis.getRow(i + 1), candidateDual[1], candidateDualDenom[1]);
1850 
1851  unsigned j = widthI[0] < widthI[1] ? 0 : 1;
1852  if (j == 0)
1853  // Subtract 1 to go from u = ceil(dual) back to floor(dual).
1854  basis.addToRow(i, i + 1, -1);
1855 
1856  // width_i(b{i+1} + u*b_i) should be minimized at our value of u.
1857  // We assert that this holds by checking that the values of width_i at
1858  // u - 1 and u + 1 are greater than or equal to the value at u. If the
1859  // width is lesser at either of the adjacent values, then our computed
1860  // value of u is clearly not the minimizer. Otherwise by convexity the
1861  // computed value of u is really the minimizer.
1862 
1863  // Check the value at u - 1.
1864  assert(gbrSimplex.computeWidth(scaleAndAddForAssert(
1865  basis.getRow(i + 1), MPInt(-1), basis.getRow(i))) >=
1866  widthI[j] &&
1867  "Computed u value does not minimize the width!");
1868  // Check the value at u + 1.
1869  assert(gbrSimplex.computeWidth(scaleAndAddForAssert(
1870  basis.getRow(i + 1), MPInt(+1), basis.getRow(i))) >=
1871  widthI[j] &&
1872  "Computed u value does not minimize the width!");
1873 
1874  dual = std::move(candidateDual[j]);
1875  dualDenom = candidateDualDenom[j];
1876  return widthI[j];
1877  }
1878 
1879  assert(i + 1 - level < width.size() && "width_{i+1} wasn't saved");
1880  // f_i(b_{i+1} + dual*b_i) == width_{i+1}(b_{i+1}) when `dual` minimizes the
1881  // LHS. (note: the basis has already been updated, so b_{i+1} + dual*b_i in
1882  // the above expression is equal to basis.getRow(i+1) below.)
1883  assert(gbrSimplex.computeWidth(basis.getRow(i + 1)) ==
1884  width[i + 1 - level]);
1885  return width[i + 1 - level];
1886  };
1887 
1888  // In the ith iteration of the loop, gbrSimplex has constraints for directions
1889  // from `level` to i - 1.
1890  unsigned i = level;
1891  while (i < basis.getNumRows() - 1) {
1892  if (i >= level + width.size()) {
1893  // We don't even know the value of f_i(b_i), so let's find that first.
1894  // We have to do this first since later we assume that width already
1895  // contains values up to and including i.
1896 
1897  assert((i == 0 || i - 1 < level + width.size()) &&
1898  "We are at level i but we don't know the value of width_{i-1}");
1899 
1900  // We don't actually use these duals at all, but it doesn't matter
1901  // because this case should only occur when i is level, and there are no
1902  // duals in that case anyway.
1903  assert(i == level && "This case should only occur when i == level");
1904  width.push_back(
1905  gbrSimplex.computeWidthAndDuals(basis.getRow(i), dual, dualDenom));
1906  }
1907 
1908  if (i >= level + dual.size()) {
1909  assert(i + 1 >= level + width.size() &&
1910  "We don't know dual_i but we know width_{i+1}");
1911  // We don't know dual for our level, so let's find it.
1912  gbrSimplex.addEqualityForDirection(basis.getRow(i));
1913  width.push_back(gbrSimplex.computeWidthAndDuals(basis.getRow(i + 1), dual,
1914  dualDenom));
1915  gbrSimplex.removeLastEquality();
1916  }
1917 
1918  // This variable stores width_i(b_{i+1} + u*b_i).
1919  Fraction widthICandidate = updateBasisWithUAndGetFCandidate(i);
1920  if (widthICandidate < epsilon * width[i - level]) {
1921  basis.swapRows(i, i + 1);
1922  width[i - level] = widthICandidate;
1923  // The values of width_{i+1}(b_{i+1}) and higher may change after the
1924  // swap, so we remove the cached values here.
1925  width.resize(i - level + 1);
1926  if (i == level) {
1927  dual.clear();
1928  continue;
1929  }
1930 
1931  gbrSimplex.removeLastEquality();
1932  i--;
1933  continue;
1934  }
1935 
1936  // Invalidate duals since the higher level needs to recompute its own duals.
1937  dual.clear();
1938  gbrSimplex.addEqualityForDirection(basis.getRow(i));
1939  i++;
1940  }
1941 }
1942 
1943 /// Search for an integer sample point using a branch and bound algorithm.
1944 ///
1945 /// Each row in the basis matrix is a vector, and the set of basis vectors
1946 /// should span the space. Initially this is the identity matrix,
1947 /// i.e., the basis vectors are just the variables.
1948 ///
1949 /// In every level, a value is assigned to the level-th basis vector, as
1950 /// follows. Compute the minimum and maximum rational values of this direction.
1951 /// If only one integer point lies in this range, constrain the variable to
1952 /// have this value and recurse to the next variable.
1953 ///
1954 /// If the range has multiple values, perform generalized basis reduction via
1955 /// reduceBasis and then compute the bounds again. Now we try constraining
1956 /// this direction in the first value in this range and "recurse" to the next
1957 /// level. If we fail to find a sample, we try assigning the direction the next
1958 /// value in this range, and so on.
1959 ///
1960 /// If no integer sample is found from any of the assignments, or if the range
1961 /// contains no integer value, then of course the polytope is empty for the
1962 /// current assignment of the values in previous levels, so we return to
1963 /// the previous level.
1964 ///
1965 /// If we reach the last level where all the variables have been assigned values
1966 /// already, then we simply return the current sample point if it is integral,
1967 /// and go back to the previous level otherwise.
1968 ///
1969 /// To avoid potentially arbitrarily large recursion depths leading to stack
1970 /// overflows, this algorithm is implemented iteratively.
1972  if (empty)
1973  return {};
1974 
1975  unsigned nDims = var.size();
1976  Matrix basis = Matrix::identity(nDims);
1977 
1978  unsigned level = 0;
1979  // The snapshot just before constraining a direction to a value at each level.
1980  SmallVector<unsigned, 8> snapshotStack;
1981  // The maximum value in the range of the direction for each level.
1982  SmallVector<MPInt, 8> upperBoundStack;
1983  // The next value to try constraining the basis vector to at each level.
1984  SmallVector<MPInt, 8> nextValueStack;
1985 
1986  snapshotStack.reserve(basis.getNumRows());
1987  upperBoundStack.reserve(basis.getNumRows());
1988  nextValueStack.reserve(basis.getNumRows());
1989  while (level != -1u) {
1990  if (level == basis.getNumRows()) {
1991  // We've assigned values to all variables. Return if we have a sample,
1992  // or go back up to the previous level otherwise.
1993  if (auto maybeSample = getSamplePointIfIntegral())
1994  return maybeSample;
1995  level--;
1996  continue;
1997  }
1998 
1999  if (level >= upperBoundStack.size()) {
2000  // We haven't populated the stack values for this level yet, so we have
2001  // just come down a level ("recursed"). Find the lower and upper bounds.
2002  // If there is more than one integer point in the range, perform
2003  // generalized basis reduction.
2004  SmallVector<MPInt, 8> basisCoeffs =
2005  llvm::to_vector<8>(basis.getRow(level));
2006  basisCoeffs.emplace_back(0);
2007 
2008  auto [minRoundedUp, maxRoundedDown] = computeIntegerBounds(basisCoeffs);
2009 
2010  // We don't have any integer values in the range.
2011  // Pop the stack and return up a level.
2012  if (minRoundedUp.isEmpty() || maxRoundedDown.isEmpty()) {
2013  assert((minRoundedUp.isEmpty() && maxRoundedDown.isEmpty()) &&
2014  "If one bound is empty, both should be.");
2015  snapshotStack.pop_back();
2016  nextValueStack.pop_back();
2017  upperBoundStack.pop_back();
2018  level--;
2019  continue;
2020  }
2021 
2022  // We already checked the empty case above.
2023  assert((minRoundedUp.isBounded() && maxRoundedDown.isBounded()) &&
2024  "Polyhedron should be bounded!");
2025 
2026  // Heuristic: if the sample point is integral at this point, just return
2027  // it.
2028  if (auto maybeSample = getSamplePointIfIntegral())
2029  return *maybeSample;
2030 
2031  if (*minRoundedUp < *maxRoundedDown) {
2032  reduceBasis(basis, level);
2033  basisCoeffs = llvm::to_vector<8>(basis.getRow(level));
2034  basisCoeffs.emplace_back(0);
2035  std::tie(minRoundedUp, maxRoundedDown) =
2036  computeIntegerBounds(basisCoeffs);
2037  }
2038 
2039  snapshotStack.push_back(getSnapshot());
2040  // The smallest value in the range is the next value to try.
2041  // The values in the optionals are guaranteed to exist since we know the
2042  // polytope is bounded.
2043  nextValueStack.push_back(*minRoundedUp);
2044  upperBoundStack.push_back(*maxRoundedDown);
2045  }
2046 
2047  assert((snapshotStack.size() - 1 == level &&
2048  nextValueStack.size() - 1 == level &&
2049  upperBoundStack.size() - 1 == level) &&
2050  "Mismatched variable stack sizes!");
2051 
2052  // Whether we "recursed" or "returned" from a lower level, we rollback
2053  // to the snapshot of the starting state at this level. (in the "recursed"
2054  // case this has no effect)
2055  rollback(snapshotStack.back());
2056  MPInt nextValue = nextValueStack.back();
2057  ++nextValueStack.back();
2058  if (nextValue > upperBoundStack.back()) {
2059  // We have exhausted the range and found no solution. Pop the stack and
2060  // return up a level.
2061  snapshotStack.pop_back();
2062  nextValueStack.pop_back();
2063  upperBoundStack.pop_back();
2064  level--;
2065  continue;
2066  }
2067 
2068  // Try the next value in the range and "recurse" into the next level.
2069  SmallVector<MPInt, 8> basisCoeffs(basis.getRow(level).begin(),
2070  basis.getRow(level).end());
2071  basisCoeffs.push_back(-nextValue);
2072  addEquality(basisCoeffs);
2073  level++;
2074  }
2075 
2076  return {};
2077 }
2078 
2079 /// Compute the minimum and maximum integer values the expression can take. We
2080 /// compute each separately.
2081 std::pair<MaybeOptimum<MPInt>, MaybeOptimum<MPInt>>
2083  MaybeOptimum<MPInt> minRoundedUp(
2085  MaybeOptimum<MPInt> maxRoundedDown(
2087  return {minRoundedUp, maxRoundedDown};
2088 }
2089 
2090 void SimplexBase::print(raw_ostream &os) const {
2091  os << "rows = " << getNumRows() << ", columns = " << getNumColumns() << "\n";
2092  if (empty)
2093  os << "Simplex marked empty!\n";
2094  os << "var: ";
2095  for (unsigned i = 0; i < var.size(); ++i) {
2096  if (i > 0)
2097  os << ", ";
2098  var[i].print(os);
2099  }
2100  os << "\ncon: ";
2101  for (unsigned i = 0; i < con.size(); ++i) {
2102  if (i > 0)
2103  os << ", ";
2104  con[i].print(os);
2105  }
2106  os << '\n';
2107  for (unsigned row = 0, e = getNumRows(); row < e; ++row) {
2108  if (row > 0)
2109  os << ", ";
2110  os << "r" << row << ": " << rowUnknown[row];
2111  }
2112  os << '\n';
2113  os << "c0: denom, c1: const";
2114  for (unsigned col = 2, e = getNumColumns(); col < e; ++col)
2115  os << ", c" << col << ": " << colUnknown[col];
2116  os << '\n';
2117  for (unsigned row = 0, numRows = getNumRows(); row < numRows; ++row) {
2118  for (unsigned col = 0, numCols = getNumColumns(); col < numCols; ++col)
2119  os << tableau(row, col) << '\t';
2120  os << '\n';
2121  }
2122  os << '\n';
2123 }
2124 
2125 void SimplexBase::dump() const { print(llvm::errs()); }
2126 
2128  if (isEmpty())
2129  return true;
2130 
2131  for (unsigned i = 0, e = rel.getNumInequalities(); i < e; ++i)
2133  return false;
2134 
2135  for (unsigned i = 0, e = rel.getNumEqualities(); i < e; ++i)
2136  if (!isRedundantEquality(rel.getEquality(i)))
2137  return false;
2138 
2139  return true;
2140 }
2141 
2142 /// Returns the type of the inequality with coefficients `coeffs`.
2143 /// Possible types are:
2144 /// Redundant The inequality is satisfied by all points in the polytope
2145 /// Cut The inequality is satisfied by some points, but not by others
2146 /// Separate The inequality is not satisfied by any point
2147 ///
2148 /// Internally, this computes the minimum and the maximum the inequality with
2149 /// coefficients `coeffs` can take. If the minimum is >= 0, the inequality holds
2150 /// for all points in the polytope, so it is redundant. If the minimum is <= 0
2151 /// and the maximum is >= 0, the points in between the minimum and the
2152 /// inequality do not satisfy it, the points in between the inequality and the
2153 /// maximum satisfy it. Hence, it is a cut inequality. If both are < 0, no
2154 /// points of the polytope satisfy the inequality, which means it is a separate
2155 /// inequality.
2158  if (minimum.isBounded() && *minimum >= Fraction(0, 1)) {
2159  return IneqType::Redundant;
2160  }
2162  if ((!minimum.isBounded() || *minimum <= Fraction(0, 1)) &&
2163  (!maximum.isBounded() || *maximum >= Fraction(0, 1))) {
2164  return IneqType::Cut;
2165  }
2166  return IneqType::Separate;
2167 }
2168 
2169 /// Checks whether the type of the inequality with coefficients `coeffs`
2170 /// is Redundant.
2172  assert(!empty &&
2173  "It is not meaningful to ask about redundancy in an empty set!");
2174  return findIneqType(coeffs) == IneqType::Redundant;
2175 }
2176 
2177 /// Check whether the equality given by `coeffs == 0` is redundant given
2178 /// the existing constraints. This is redundant when `coeffs` is already
2179 /// always zero under the existing constraints. `coeffs` is always zero
2180 /// when the minimum and maximum value that `coeffs` can take are both zero.
2182  assert(!empty &&
2183  "It is not meaningful to ask about redundancy in an empty set!");
2186  assert((!minimum.isEmpty() && !maximum.isEmpty()) &&
2187  "Optima should be non-empty for a non-empty set");
2188  return minimum.isBounded() && maximum.isBounded() &&
2189  *maximum == Fraction(0, 1) && *minimum == Fraction(0, 1);
2190 }
static Value max(ImplicitLocOpBuilder &builder, Value value, Value bound)
static LLVM_ATTRIBUTE_UNUSED SmallVector< MPInt, 8 > scaleAndAddForAssert(ArrayRef< MPInt > a, const MPInt &scale, ArrayRef< MPInt > b)
Definition: Simplex.cpp:26
static bool isRangeDivisibleBy(ArrayRef< MPInt > range, const MPInt &divisor)
Definition: Simplex.cpp:354
const int nullIndex
Definition: Simplex.cpp:21
An IntegerRelation represents the set of points from a PresburgerSpace that satisfy a list of affine ...
ArrayRef< MPInt > getEquality(unsigned idx) const
void addInequality(ArrayRef< MPInt > inEq)
Adds an inequality (>= 0) from the coefficients specified in inEq.
void truncate(const CountsSnapshot &counts)
void addLocalFloorDiv(ArrayRef< MPInt > dividend, const MPInt &divisor)
Adds a new local variable as the floordiv of an affine function of other variables,...
ArrayRef< MPInt > getInequality(unsigned idx) const
DivisionRepr getLocalReprs(std::vector< MaybeLocalRepr > *repr=nullptr) const
Returns a DivisonRepr representing the division representation of local variables in the constraint s...
void undoLastConstraint() final
Undo the addition of the last constraint.
Definition: Simplex.cpp:1191
LogicalResult moveRowUnknownToColumn(unsigned row)
Try to move the specified row to column orientation while preserving the lexicopositivity of the basi...
Definition: Simplex.cpp:758
void addInequality(ArrayRef< MPInt > coeffs) final
Add an inequality to the tableau.
Definition: Simplex.cpp:1569
LogicalResult addCut(unsigned row)
Given a row that has a non-integer sample value, add an inequality to cut away this fractional sample...
Definition: Simplex.cpp:263
unsigned getLexMinPivotColumn(unsigned row, unsigned colA, unsigned colB) const
Given two potential pivot columns for a row, return the one that results in the lexicographically sma...
Definition: Simplex.cpp:774
unsigned getSnapshot()
Get a snapshot of the current state. This is used for rolling back.
Definition: Simplex.h:429
void appendSymbol()
Add new symbolic variables to the end of the list of variables.
Definition: Simplex.cpp:347
MaybeOptimum< SmallVector< MPInt, 8 > > findIntegerLexMin()
Return the lexicographically minimum integer solution to the constraints.
Definition: Simplex.cpp:288
bool isRedundantInequality(ArrayRef< MPInt > coeffs)
Definition: Simplex.cpp:325
MaybeOptimum< SmallVector< Fraction, 8 > > findRationalLexMin()
Return the lexicographically minimum rational solution to the constraints.
Definition: Simplex.cpp:217
bool isSeparateInequality(ArrayRef< MPInt > coeffs)
Return whether the specified inequality is redundant/separate for the polytope.
Definition: Simplex.cpp:319
This class provides support for multi-precision arithmetic.
Definition: MPInt.h:87
This is a class to represent a resizable matrix.
Definition: Matrix.h:35
void swapColumns(unsigned column, unsigned otherColumn)
Swap the given columns.
Definition: Matrix.cpp:78
unsigned appendExtraRow()
Add an extra row at the bottom of the matrix and return its position.
Definition: Matrix.cpp:39
static Matrix identity(unsigned dimension)
Return the identity matrix of the specified dimension.
Definition: Matrix.cpp:24
MPInt normalizeRow(unsigned row, unsigned nCols)
Divide the first nCols of the specified row by their GCD.
Definition: Matrix.cpp:225
void resizeVertically(unsigned newNRows)
Definition: Matrix.cpp:64
unsigned getNumRows() const
Definition: Matrix.h:76
void swapRows(unsigned row, unsigned otherRow)
Swap the given rows.
Definition: Matrix.cpp:69
void resizeHorizontally(unsigned newNColumns)
Definition: Matrix.cpp:52
void addToRow(unsigned sourceRow, unsigned targetRow, const MPInt &scale)
Add scale multiples of the source row to the target row.
Definition: Matrix.cpp:194
void reserveRows(unsigned rows)
Reserve enough space to resize to the specified number of rows without reallocations.
Definition: Matrix.cpp:35
MutableArrayRef< MPInt > getRow(unsigned row)
Get a [Mutable]ArrayRef corresponding to the specified row.
Definition: Matrix.cpp:87
bool isBounded() const
Definition: Utils.h:51
bool isUnbounded() const
Definition: Utils.h:52
This class represents a multi-affine function with the domain as Z^d, where d is the number of domain...
Definition: PWMAFunction.h:36
const PresburgerSpace & getSpace() const
Definition: PWMAFunction.h:152
void addPiece(const Piece &piece)
unsigned getNumOutputs() const
Definition: PWMAFunction.h:162
void unionInPlace(const IntegerRelation &disjunct)
Mutate this set, turning it into the union of this set and the given disjunct.
PresburgerSpace is the space of all possible values of a tuple of integer valued variables/variables.
static PresburgerSpace getRelationSpace(unsigned numDomain=0, unsigned numRange=0, unsigned numSymbols=0, unsigned numLocals=0)
unsigned insertVar(VarKind kind, unsigned pos, unsigned num=1)
Insert num variables of the specified kind at position pos.
The Simplex class implements a version of the Simplex and Generalized Basis Reduction algorithms,...
Definition: Simplex.h:157
Optional< unsigned > findAnyPivotRow(unsigned col)
Return any row that this column can be pivoted with, ignoring tableau consistency.
Definition: Simplex.cpp:1154
unsigned addZeroRow(bool makeRestricted=false)
Add a new row to the tableau and the associated data structures.
Definition: Simplex.cpp:90
bool isEmpty() const
Returns true if the tableau is empty (has conflicting constraints), false otherwise.
Definition: Simplex.cpp:1051
void appendVariable(unsigned count=1)
Add new variables to the end of the list of variables.
Definition: Simplex.cpp:1299
virtual void undoLastConstraint()=0
Undo the addition of the last constraint.
SmallVector< int, 8 > rowUnknown
These hold the indexes of the unknown at a given row or column position.
Definition: Simplex.h:362
SmallVector< SmallVector< int, 8 >, 8 > savedBases
Holds a vector of bases.
Definition: Simplex.h:353
void intersectIntegerRelation(const IntegerRelation &rel)
Add all the constraints from the given IntegerRelation.
Definition: Simplex.cpp:1314
SmallVector< UndoLogEntry, 8 > undoLog
Holds a log of operations, used for rolling back to a previous state.
Definition: Simplex.h:348
bool usingBigM
Stores whether or not a big M column is present in the tableau.
Definition: Simplex.h:330
unsigned getSnapshot() const
Get a snapshot of the current state.
Definition: Simplex.cpp:1117
void print(raw_ostream &os) const
Print the tableau's internal state.
Definition: Simplex.cpp:2090
UndoLogEntry
Enum to denote operations that need to be undone during rollback.
Definition: Simplex.h:305
unsigned getNumRows() const
Definition: Simplex.h:326
unsigned addRow(ArrayRef< MPInt > coeffs, bool makeRestricted=false)
Add a new row to the tableau and the associated data structures.
Definition: Simplex.cpp:104
const Unknown & unknownFromRow(unsigned row) const
Returns the unknown associated with row.
Definition: Simplex.cpp:70
SmallVector< int, 8 > colUnknown
Definition: Simplex.h:362
SmallVector< Unknown, 8 > var
Definition: Simplex.h:365
virtual void addInequality(ArrayRef< MPInt > coeffs)=0
Add an inequality to the tableau.
unsigned getSnapshotBasis()
Get a snapshot of the current state including the basis.
Definition: Simplex.cpp:1119
Matrix tableau
The matrix representing the tableau.
Definition: Simplex.h:341
unsigned getNumFixedCols() const
Return the number of fixed columns, as described in the constructor above, this is the number of colu...
Definition: Simplex.h:325
SmallVector< Unknown, 8 > con
These hold information about each unknown.
Definition: Simplex.h:365
void markEmpty()
Mark the tableau as being empty.
Definition: Simplex.cpp:1074
bool empty
This is true if the tableau has been detected to be empty, false otherwise.
Definition: Simplex.h:345
void swapColumns(unsigned i, unsigned j)
Definition: Simplex.cpp:1062
void removeLastConstraintRowOrientation()
Remove the last constraint, which must be in row orientation.
Definition: Simplex.cpp:1131
void addEquality(ArrayRef< MPInt > coeffs)
Add an equality to the tableau.
Definition: Simplex.cpp:1104
const Unknown & unknownFromColumn(unsigned col) const
Returns the unknown associated with col.
Definition: Simplex.cpp:65
void rollback(unsigned snapshot)
Rollback to a snapshot. This invalidates all later snapshots.
Definition: Simplex.cpp:1269
void pivot(unsigned row, unsigned col)
Pivot the row with the column.
Definition: Simplex.cpp:932
void swapRows(unsigned i, unsigned j)
Swap the two rows/columns in the tableau and associated data structures.
Definition: Simplex.cpp:1053
void undo(UndoLogEntry entry)
Undo the operation represented by the log entry.
Definition: Simplex.cpp:1207
void addDivisionVariable(ArrayRef< MPInt > coeffs, const MPInt &denom)
Append a new variable to the simplex and constrain it such that its only integer value is the floor d...
Definition: Simplex.cpp:1282
const Unknown & unknownFromIndex(int index) const
Returns the unknown associated with index.
Definition: Simplex.cpp:60
unsigned nSymbol
The number of parameters.
Definition: Simplex.h:338
unsigned nRedundant
The number of redundant rows in the tableau.
Definition: Simplex.h:334
unsigned getNumVariables() const
Returns the number of variables in the tableau.
Definition: Simplex.cpp:1112
void swapRowWithCol(unsigned row, unsigned col)
Swap the row with the column in the tableau's data structures but not the tableau itself.
Definition: Simplex.cpp:895
unsigned getNumColumns() const
Definition: Simplex.h:327
unsigned getNumConstraints() const
Returns the number of constraints in the tableau.
Definition: Simplex.cpp:1113
Takes a snapshot of the simplex state on construction and rolls back to the snapshot on destruction.
Definition: Simplex.h:872
The Simplex class uses the Normal pivot rule and supports integer emptiness checks as well as detecti...
Definition: Simplex.h:695
MaybeOptimum< Fraction > computeOptimum(Direction direction, ArrayRef< MPInt > coeffs)
Compute the maximum or minimum value of the given expression, depending on direction.
Definition: Simplex.cpp:1342
bool isMarkedRedundant(unsigned constraintIndex) const
Returns whether the specified constraint has been marked as redundant.
Definition: Simplex.cpp:1387
void addInequality(ArrayRef< MPInt > coeffs) final
Add an inequality to the tableau.
Definition: Simplex.cpp:1091
Optional< SmallVector< MPInt, 8 > > getSamplePointIfIntegral() const
Returns the current sample point if it is integral.
Definition: Simplex.cpp:1603
bool isRedundantEquality(ArrayRef< MPInt > coeffs)
Check if the specified equality already holds in the polytope.
Definition: Simplex.cpp:2181
static Simplex makeProduct(const Simplex &a, const Simplex &b)
Make a tableau to represent a pair of points in the given tableaus, one in tableau A and one in B.
Definition: Simplex.cpp:1476
MaybeOptimum< Fraction > computeRowOptimum(Direction direction, unsigned row)
Compute the maximum or minimum value of the given row, depending on direction.
Definition: Simplex.cpp:1323
bool isRationalSubsetOf(const IntegerRelation &rel)
Returns true if this Simplex's polytope is a rational subset of rel.
Definition: Simplex.cpp:2127
std::pair< MaybeOptimum< MPInt >, MaybeOptimum< MPInt > > computeIntegerBounds(ArrayRef< MPInt > coeffs)
Returns a (min, max) pair denoting the minimum and maximum integer values of the given expression.
Definition: Simplex.cpp:2082
bool isRedundantInequality(ArrayRef< MPInt > coeffs)
Check if the specified inequality already holds in the polytope.
Definition: Simplex.cpp:2171
bool isBoundedAlongConstraint(unsigned constraintIndex)
Returns whether the perpendicular of the specified constraint is a is a direction along which the pol...
Definition: Simplex.cpp:1377
bool isUnbounded()
Returns true if the polytope is unbounded, i.e., extends to infinity in some direction.
Definition: Simplex.cpp:1448
IneqType findIneqType(ArrayRef< MPInt > coeffs)
Returns the type of the inequality with coefficients coeffs.
Definition: Simplex.cpp:2156
Optional< SmallVector< Fraction, 8 > > getRationalSample() const
Returns the current sample point, which may contain non-integer (rational) coordinates.
Definition: Simplex.cpp:1548
Optional< SmallVector< MPInt, 8 > > findIntegerSample()
Returns an integer sample point if one exists, or std::nullopt otherwise.
Definition: Simplex.cpp:1971
SymbolicLexMin computeSymbolicIntegerLexMin()
The lexmin will be stored as a function lexmin from symbols to non-symbols in the result.
Definition: Simplex.cpp:518
Given a simplex for a polytope, construct a new simplex whose variables are identified with a pair of...
Definition: Simplex.cpp:1631
Fraction computeWidth(ArrayRef< MPInt > dir)
Compute max(dotProduct(dir, x - y)).
Definition: Simplex.cpp:1649
Fraction computeWidthAndDuals(ArrayRef< MPInt > dir, SmallVectorImpl< MPInt > &dual, MPInt &dualDenom)
Compute max(dotProduct(dir, x - y)) and save the dual variables for only the direction equalities to ...
Definition: Simplex.cpp:1658
void removeLastEquality()
Remove the last equality that was added through addEqualityForDirection.
Definition: Simplex.cpp:1717
void addEqualityForDirection(ArrayRef< MPInt > dir)
Add an equality dotProduct(dir, x - y) == 0.
Definition: Simplex.cpp:1642
GBRSimplex(const Simplex &originalSimplex)
Definition: Simplex.cpp:1635
SmallVector< AffineExpr, 4 > concat(ArrayRef< AffineExpr > a, ArrayRef< AffineExpr > b)
Return the vector that is the concatenation of a and b.
Definition: LinalgOps.cpp:1679
LLVM_ATTRIBUTE_ALWAYS_INLINE MPInt mod(const MPInt &lhs, const MPInt &rhs)
is always non-negative.
Definition: MPInt.h:393
void normalizeDiv(MutableArrayRef< MPInt > num, MPInt &denom)
Normalize the given (numerator, denominator) pair by dividing out the common factors between them.
Definition: Utils.cpp:355
MPInt ceil(const Fraction &f)
Definition: Fraction.h:70
MPInt normalizeRange(MutableArrayRef< MPInt > range)
Divide the range by its gcd and return the gcd.
Definition: Utils.cpp:346
MPInt floor(const Fraction &f)
Definition: Fraction.h:68
SmallVector< MPInt, 8 > getComplementIneq(ArrayRef< MPInt > ineq)
Return the complement of the given inequality.
Definition: Utils.cpp:371
LLVM_ATTRIBUTE_ALWAYS_INLINE MPInt lcm(const MPInt &a, const MPInt &b)
Returns the least common multiple of 'a' and 'b'.
Definition: MPInt.h:407
LLVM_ATTRIBUTE_ALWAYS_INLINE MPInt floorDiv(const MPInt &lhs, const MPInt &rhs)
Definition: MPInt.h:382
Include the generated interface declarations.
LogicalResult failure(bool isFailure=true)
Utility function to generate a LogicalResult.
Definition: LogicalResult.h:62
LogicalResult success(bool isSuccess=true)
Utility function to generate a LogicalResult.
Definition: LogicalResult.h:56
bool failed(LogicalResult result)
Utility function that returns true if the provided LogicalResult corresponds to a failure value.
Definition: LogicalResult.h:72
This class represents an efficient way to signal success or failure.
Definition: LogicalResult.h:26
bool failed() const
Returns true if the provided LogicalResult corresponds to a failure value.
Definition: LogicalResult.h:44
A class to represent fractions.
Definition: Fraction.h:28
MPInt getAsInteger() const
Definition: Fraction.h:46
The struct CountsSnapshot stores the count of each VarKind, and also of each constraint type.
An Unknown is either a variable or a constraint.
Definition: Simplex.h:238
Represents the result of a symbolic lexicographic minimization computation.
Definition: Simplex.h:533
PWMAFunction lexmin
This maps assignments of symbols to the corresponding lexmin.
Definition: Simplex.h:541
PresburgerSet unboundedDomain
Contains all assignments to the symbols that made the lexmin unbounded.
Definition: Simplex.h:545
Eliminates variable at the specified position using Fourier-Motzkin variable elimination.