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