//===- Simplex.cpp - MLIR Simplex Class -----------------------------------===// // // Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions. // See https://llvm.org/LICENSE.txt for license information. // SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception // //===----------------------------------------------------------------------===// #include "mlir/Analysis/Presburger/Simplex.h" #include "mlir/Analysis/Presburger/Fraction.h" #include "mlir/Analysis/Presburger/IntegerRelation.h" #include "mlir/Analysis/Presburger/Matrix.h" #include "mlir/Analysis/Presburger/PresburgerSpace.h" #include "mlir/Analysis/Presburger/Utils.h" #include "llvm/ADT/DynamicAPInt.h" #include "llvm/ADT/STLExtras.h" #include "llvm/ADT/SmallBitVector.h" #include "llvm/ADT/SmallVector.h" #include "llvm/Support/Compiler.h" #include "llvm/Support/ErrorHandling.h" #include "llvm/Support/LogicalResult.h" #include "llvm/Support/raw_ostream.h" #include <cassert> #include <functional> #include <limits> #include <optional> #include <tuple> #include <utility> usingnamespacemlir; usingnamespacepresburger; Direction; const int nullIndex = …; // Return a + scale*b; LLVM_ATTRIBUTE_UNUSED static SmallVector<DynamicAPInt, 8> scaleAndAddForAssert(ArrayRef<DynamicAPInt> a, const DynamicAPInt &scale, ArrayRef<DynamicAPInt> b) { … } SimplexBase::SimplexBase(unsigned nVar, bool mustUseBigM) : … { … } SimplexBase::SimplexBase(unsigned nVar, bool mustUseBigM, const llvm::SmallBitVector &isSymbol) : … { … } const Simplex::Unknown &SimplexBase::unknownFromIndex(int index) const { … } const Simplex::Unknown &SimplexBase::unknownFromColumn(unsigned col) const { … } const Simplex::Unknown &SimplexBase::unknownFromRow(unsigned row) const { … } Simplex::Unknown &SimplexBase::unknownFromIndex(int index) { … } Simplex::Unknown &SimplexBase::unknownFromColumn(unsigned col) { … } Simplex::Unknown &SimplexBase::unknownFromRow(unsigned row) { … } unsigned SimplexBase::addZeroRow(bool makeRestricted) { … } /// Add a new row to the tableau corresponding to the given constant term and /// list of coefficients. The coefficients are specified as a vector of /// (variable index, coefficient) pairs. unsigned SimplexBase::addRow(ArrayRef<DynamicAPInt> coeffs, bool makeRestricted) { … } namespace { bool signMatchesDirection(const DynamicAPInt &elem, Direction direction) { … } Direction flippedDirection(Direction direction) { … } } // namespace /// We simply make the tableau consistent while maintaining a lexicopositive /// basis transform, and then return the sample value. If the tableau becomes /// empty, we return empty. /// /// Let the variables be x = (x_1, ... x_n). /// Let the basis unknowns be y = (y_1, ... y_n). /// We have that x = A*y + b for some n x n matrix A and n x 1 column vector b. /// /// As we will show below, A*y is either zero or lexicopositive. /// Adding a lexicopositive vector to b will make it lexicographically /// greater, so A*y + b is always equal to or lexicographically greater than b. /// Thus, since we can attain x = b, that is the lexicographic minimum. /// /// We have that every column in A is lexicopositive, i.e., has at least /// one non-zero element, with the first such element being positive. Since for /// the tableau to be consistent we must have non-negative sample values not /// only for the constraints but also for the variables, we also have x >= 0 and /// y >= 0, by which we mean every element in these vectors is non-negative. /// /// Proof that if every column in A is lexicopositive, and y >= 0, then /// A*y is zero or lexicopositive. Begin by considering A_1, the first row of A. /// If this row is all zeros, then (A*y)_1 = (A_1)*y = 0; proceed to the next /// row. If we run out of rows, A*y is zero and we are done; otherwise, we /// encounter some row A_i that has a non-zero element. Every column is /// lexicopositive and so has some positive element before any negative elements /// occur, so the element in this row for any column, if non-zero, must be /// positive. Consider (A*y)_i = (A_i)*y. All the elements in both vectors are /// non-negative, so if this is non-zero then it must be positive. Then the /// first non-zero element of A*y is positive so A*y is lexicopositive. /// /// Otherwise, if (A_i)*y is zero, then for every column j that had a non-zero /// element in A_i, y_j is zero. Thus these columns have no contribution to A*y /// and we can completely ignore these columns of A. We now continue downwards, /// looking for rows of A that have a non-zero element other than in the ignored /// columns. If we find one, say A_k, once again these elements must be positive /// since they are the first non-zero element in each of these columns, so if /// (A_k)*y is not zero then we have that A*y is lexicopositive and if not we /// add these to the set of ignored columns and continue to the next row. If we /// run out of rows, then A*y is zero and we are done. MaybeOptimum<SmallVector<Fraction, 8>> LexSimplex::findRationalLexMin() { … } /// Given a row that has a non-integer sample value, add an inequality such /// that this fractional sample value is cut away from the polytope. The added /// inequality will be such that no integer points are removed. i.e., the /// integer lexmin, if it exists, is the same with and without this constraint. /// /// Let the row be /// (c + coeffM*M + a_1*s_1 + ... + a_m*s_m + b_1*y_1 + ... + b_n*y_n)/d, /// where s_1, ... s_m are the symbols and /// y_1, ... y_n are the other basis unknowns. /// /// For this to be an integer, we want /// coeffM*M + a_1*s_1 + ... + a_m*s_m + b_1*y_1 + ... + b_n*y_n = -c (mod d) /// Note that this constraint must always hold, independent of the basis, /// becuse the row unknown's value always equals this expression, even if *we* /// later compute the sample value from a different expression based on a /// different basis. /// /// Let us assume that M has a factor of d in it. Imposing this constraint on M /// does not in any way hinder us from finding a value of M that is big enough. /// Moreover, this function is only called when the symbolic part of the sample, /// a_1*s_1 + ... + a_m*s_m, is known to be an integer. /// /// Also, we can safely reduce the coefficients modulo d, so we have: /// /// (b_1%d)y_1 + ... + (b_n%d)y_n = (-c%d) + k*d for some integer `k` /// /// Note that all coefficient modulos here are non-negative. Also, all the /// unknowns are non-negative here as both constraints and variables are /// non-negative in LexSimplexBase. (We used the big M trick to make the /// variables non-negative). Therefore, the LHS here is non-negative. /// Since 0 <= (-c%d) < d, k is the quotient of dividing the LHS by d and /// is therefore non-negative as well. /// /// So we have /// ((b_1%d)y_1 + ... + (b_n%d)y_n - (-c%d))/d >= 0. /// /// The constraint is violated when added (it would be useless otherwise) /// so we immediately try to move it to a column. LogicalResult LexSimplexBase::addCut(unsigned row) { … } std::optional<unsigned> LexSimplex::maybeGetNonIntegralVarRow() const { … } MaybeOptimum<SmallVector<DynamicAPInt, 8>> LexSimplex::findIntegerLexMin() { … } bool LexSimplex::isSeparateInequality(ArrayRef<DynamicAPInt> coeffs) { … } bool LexSimplex::isRedundantInequality(ArrayRef<DynamicAPInt> coeffs) { … } SmallVector<DynamicAPInt, 8> SymbolicLexSimplex::getSymbolicSampleNumerator(unsigned row) const { … } SmallVector<DynamicAPInt, 8> SymbolicLexSimplex::getSymbolicSampleIneq(unsigned row) const { … } void LexSimplexBase::appendSymbol() { … } static bool isRangeDivisibleBy(ArrayRef<DynamicAPInt> range, const DynamicAPInt &divisor) { … } bool SymbolicLexSimplex::isSymbolicSampleIntegral(unsigned row) const { … } /// This proceeds similarly to LexSimplexBase::addCut(). We are given a row that /// has a symbolic sample value with fractional coefficients. /// /// Let the row be /// (c + coeffM*M + sum_i a_i*s_i + sum_j b_j*y_j)/d, /// where s_1, ... s_m are the symbols and /// y_1, ... y_n are the other basis unknowns. /// /// As in LexSimplex::addCut, for this to be an integer, we want /// /// coeffM*M + sum_j b_j*y_j = -c + sum_i (-a_i*s_i) (mod d) /// /// This time, a_1*s_1 + ... + a_m*s_m may not be an integer. We find that /// /// sum_i (b_i%d)y_i = ((-c%d) + sum_i (-a_i%d)s_i)%d + k*d for some integer k /// /// where we take a modulo of the whole symbolic expression on the right to /// bring it into the range [0, d - 1]. Therefore, as in addCut(), /// k is the quotient on dividing the LHS by d, and since LHS >= 0, we have /// k >= 0 as well. If all the a_i are divisible by d, then we can add the /// constraint directly. Otherwise, we realize the modulo of the symbolic /// expression by adding a division variable /// /// q = ((-c%d) + sum_i (-a_i%d)s_i)/d /// /// to the symbol domain, so the equality becomes /// /// sum_i (b_i%d)y_i = (-c%d) + sum_i (-a_i%d)s_i - q*d + k*d for some integer k /// /// So the cut is /// (sum_i (b_i%d)y_i - (-c%d) - sum_i (-a_i%d)s_i + q*d)/d >= 0 /// This constraint is violated when added so we immediately try to move it to a /// column. LogicalResult SymbolicLexSimplex::addSymbolicCut(unsigned row) { … } void SymbolicLexSimplex::recordOutput(SymbolicLexOpt &result) const { … } std::optional<unsigned> SymbolicLexSimplex::maybeGetAlwaysViolatedRow() { … } std::optional<unsigned> SymbolicLexSimplex::maybeGetNonIntegralVarRow() { … } /// The non-branching pivots are just the ones moving the rows /// that are always violated in the symbol domain. LogicalResult SymbolicLexSimplex::doNonBranchingPivots() { … } SymbolicLexOpt SymbolicLexSimplex::computeSymbolicIntegerLexMin() { … } bool LexSimplex::rowIsViolated(unsigned row) const { … } std::optional<unsigned> LexSimplex::maybeGetViolatedRow() const { … } /// We simply look for violated rows and keep trying to move them to column /// orientation, which always succeeds unless the constraints have no solution /// in which case we just give up and return. LogicalResult LexSimplex::restoreRationalConsistency() { … } // Move the row unknown to column orientation while preserving lexicopositivity // of the basis transform. The sample value of the row must be non-positive. // // We only consider pivots where the pivot element is positive. Suppose no such // pivot exists, i.e., some violated row has no positive coefficient for any // basis unknown. The row can be represented as (s + c_1*u_1 + ... + c_n*u_n)/d, // where d is the denominator, s is the sample value and the c_i are the basis // coefficients. If s != 0, then since any feasible assignment of the basis // satisfies u_i >= 0 for all i, and we have s < 0 as well as c_i < 0 for all i, // any feasible assignment would violate this row and therefore the constraints // have no solution. // // We can preserve lexicopositivity by picking the pivot column with positive // pivot element that makes the lexicographically smallest change to the sample // point. // // Proof. Let // x = (x_1, ... x_n) be the variables, // z = (z_1, ... z_m) be the constraints, // y = (y_1, ... y_n) be the current basis, and // define w = (x_1, ... x_n, z_1, ... z_m) = B*y + s. // B is basically the simplex tableau of our implementation except that instead // of only describing the transform to get back the non-basis unknowns, it // defines the values of all the unknowns in terms of the basis unknowns. // Similarly, s is the column for the sample value. // // Our goal is to show that each column in B, restricted to the first n // rows, is lexicopositive after the pivot if it is so before. This is // equivalent to saying the columns in the whole matrix are lexicopositive; // there must be some non-zero element in every column in the first n rows since // the n variables cannot be spanned without using all the n basis unknowns. // // Consider a pivot where z_i replaces y_j in the basis. Recall the pivot // transform for the tableau derived for SimplexBase::pivot: // // pivot col other col pivot col other col // pivot row a b -> pivot row 1/a -b/a // other row c d other row c/a d - bc/a // // Similarly, a pivot results in B changing to B' and c to c'; the difference // between the tableau and these matrices B and B' is that there is no special // case for the pivot row, since it continues to represent the same unknown. The // same formula applies for all rows: // // B'.col(j) = B.col(j) / B(i,j) // B'.col(k) = B.col(k) - B(i,k) * B.col(j) / B(i,j) for k != j // and similarly, s' = s - s_i * B.col(j) / B(i,j). // // If s_i == 0, then the sample value remains unchanged. Otherwise, if s_i < 0, // the change in sample value when pivoting with column a is lexicographically // smaller than that when pivoting with column b iff B.col(a) / B(i, a) is // lexicographically smaller than B.col(b) / B(i, b). // // Since B(i, j) > 0, column j remains lexicopositive. // // For the other columns, suppose C.col(k) is not lexicopositive. // This means that for some p, for all t < p, // C(t,k) = 0 => B(t,k) = B(t,j) * B(i,k) / B(i,j) and // C(t,k) < 0 => B(p,k) < B(t,j) * B(i,k) / B(i,j), // which is in contradiction to the fact that B.col(j) / B(i,j) must be // lexicographically smaller than B.col(k) / B(i,k), since it lexicographically // minimizes the change in sample value. LogicalResult LexSimplexBase::moveRowUnknownToColumn(unsigned row) { … } unsigned LexSimplexBase::getLexMinPivotColumn(unsigned row, unsigned colA, unsigned colB) const { … } /// Find a pivot to change the sample value of the row in the specified /// direction. The returned pivot row will involve `row` if and only if the /// unknown is unbounded in the specified direction. /// /// To increase (resp. decrease) the value of a row, we need to find a live /// column with a non-zero coefficient. If the coefficient is positive, we need /// to increase (decrease) the value of the column, and if the coefficient is /// negative, we need to decrease (increase) the value of the column. Also, /// we cannot decrease the sample value of restricted columns. /// /// If multiple columns are valid, we break ties by considering a lexicographic /// ordering where we prefer unknowns with lower index. std::optional<SimplexBase::Pivot> Simplex::findPivot(int row, Direction direction) const { … } /// Swap the associated unknowns for the row and the column. /// /// First we swap the index associated with the row and column. Then we update /// the unknowns to reflect their new position and orientation. void SimplexBase::swapRowWithCol(unsigned row, unsigned col) { … } void SimplexBase::pivot(Pivot pair) { … } /// Pivot pivotRow and pivotCol. /// /// Let R be the pivot row unknown and let C be the pivot col unknown. /// Since initially R = a*C + sum b_i * X_i /// (where the sum is over the other column's unknowns, x_i) /// C = (R - (sum b_i * X_i))/a /// /// Let u be some other row unknown. /// u = c*C + sum d_i * X_i /// So u = c*(R - sum b_i * X_i)/a + sum d_i * X_i /// /// This results in the following transform: /// pivot col other col pivot col other col /// pivot row a b -> pivot row 1/a -b/a /// other row c d other row c/a d - bc/a /// /// Taking into account the common denominators p and q: /// /// pivot col other col pivot col other col /// pivot row a/p b/p -> pivot row p/a -b/a /// other row c/q d/q other row cp/aq (da - bc)/aq /// /// The pivot row transform is accomplished be swapping a with the pivot row's /// common denominator and negating the pivot row except for the pivot column /// element. void SimplexBase::pivot(unsigned pivotRow, unsigned pivotCol) { … } /// Perform pivots until the unknown has a non-negative sample value or until /// no more upward pivots can be performed. Return success if we were able to /// bring the row to a non-negative sample value, and failure otherwise. LogicalResult Simplex::restoreRow(Unknown &u) { … } /// Find a row that can be used to pivot the column in the specified direction. /// This returns an empty optional if and only if the column is unbounded in the /// specified direction (ignoring skipRow, if skipRow is set). /// /// If skipRow is set, this row is not considered, and (if it is restricted) its /// restriction may be violated by the returned pivot. Usually, skipRow is set /// because we don't want to move it to column position unless it is unbounded, /// and we are either trying to increase the value of skipRow or explicitly /// trying to make skipRow negative, so we are not concerned about this. /// /// If the direction is up (resp. down) and a restricted row has a negative /// (positive) coefficient for the column, then this row imposes a bound on how /// much the sample value of the column can change. Such a row with constant /// term c and coefficient f for the column imposes a bound of c/|f| on the /// change in sample value (in the specified direction). (note that c is /// non-negative here since the row is restricted and the tableau is consistent) /// /// We iterate through the rows and pick the row which imposes the most /// stringent bound, since pivoting with a row changes the row's sample value to /// 0 and hence saturates the bound it imposes. We break ties between rows that /// impose the same bound by considering a lexicographic ordering where we /// prefer unknowns with lower index value. std::optional<unsigned> Simplex::findPivotRow(std::optional<unsigned> skipRow, Direction direction, unsigned col) const { … } bool SimplexBase::isEmpty() const { … } void SimplexBase::swapRows(unsigned i, unsigned j) { … } void SimplexBase::swapColumns(unsigned i, unsigned j) { … } /// Mark this tableau empty and push an entry to the undo stack. void SimplexBase::markEmpty() { … } /// Add an inequality to the tableau. If coeffs is c_0, c_1, ... c_n, where n /// is the current number of variables, then the corresponding inequality is /// c_n + c_0*x_0 + c_1*x_1 + ... + c_{n-1}*x_{n-1} >= 0. /// /// We add the inequality and mark it as restricted. We then try to make its /// sample value non-negative. If this is not possible, the tableau has become /// empty and we mark it as such. void Simplex::addInequality(ArrayRef<DynamicAPInt> coeffs) { … } /// Add an equality to the tableau. If coeffs is c_0, c_1, ... c_n, where n /// is the current number of variables, then the corresponding equality is /// c_n + c_0*x_0 + c_1*x_1 + ... + c_{n-1}*x_{n-1} == 0. /// /// We simply add two opposing inequalities, which force the expression to /// be zero. void SimplexBase::addEquality(ArrayRef<DynamicAPInt> coeffs) { … } unsigned SimplexBase::getNumVariables() const { … } unsigned SimplexBase::getNumConstraints() const { … } /// Return a snapshot of the current state. This is just the current size of the /// undo log. unsigned SimplexBase::getSnapshot() const { … } unsigned SimplexBase::getSnapshotBasis() { … } void SimplexBase::removeLastConstraintRowOrientation() { … } // This doesn't find a pivot row only if the column has zero // coefficients for every row. // // If the unknown is a constraint, this can't happen, since it was added // initially as a row. Such a row could never have been pivoted to a column. So // a pivot row will always be found if we have a constraint. // // If we have a variable, then the column has zero coefficients for every row // iff no constraints have been added with a non-zero coefficient for this row. std::optional<unsigned> SimplexBase::findAnyPivotRow(unsigned col) { … } // It's not valid to remove the constraint by deleting the column since this // would result in an invalid basis. void Simplex::undoLastConstraint() { … } // It's not valid to remove the constraint by deleting the column since this // would result in an invalid basis. void LexSimplexBase::undoLastConstraint() { … } void SimplexBase::undo(UndoLogEntry entry) { … } /// Rollback to the specified snapshot. /// /// We undo all the log entries until the log size when the snapshot was taken /// is reached. void SimplexBase::rollback(unsigned snapshot) { … } /// We add the usual floor division constraints: /// `0 <= coeffs - denom*q <= denom - 1`, where `q` is the new division /// variable. /// /// This constrains the remainder `coeffs - denom*q` to be in the /// range `[0, denom - 1]`, which fixes the integer value of the quotient `q`. void SimplexBase::addDivisionVariable(ArrayRef<DynamicAPInt> coeffs, const DynamicAPInt &denom) { … } void SimplexBase::appendVariable(unsigned count) { … } /// Add all the constraints from the given IntegerRelation. void SimplexBase::intersectIntegerRelation(const IntegerRelation &rel) { … } MaybeOptimum<Fraction> Simplex::computeRowOptimum(Direction direction, unsigned row) { … } /// Compute the optimum of the specified expression in the specified direction, /// or std::nullopt if it is unbounded. MaybeOptimum<Fraction> Simplex::computeOptimum(Direction direction, ArrayRef<DynamicAPInt> coeffs) { … } MaybeOptimum<Fraction> Simplex::computeOptimum(Direction direction, Unknown &u) { … } bool Simplex::isBoundedAlongConstraint(unsigned constraintIndex) { … } /// Redundant constraints are those that are in row orientation and lie in /// rows 0 to nRedundant - 1. bool Simplex::isMarkedRedundant(unsigned constraintIndex) const { … } /// Mark the specified row redundant. /// /// This is done by moving the unknown to the end of the block of redundant /// rows (namely, to row nRedundant) and incrementing nRedundant to /// accomodate the new redundant row. void Simplex::markRowRedundant(Unknown &u) { … } /// Find a subset of constraints that is redundant and mark them redundant. void Simplex::detectRedundant(unsigned offset, unsigned count) { … } bool Simplex::isUnbounded() { … } /// Make a tableau to represent a pair of points in the original tableau. /// /// The product constraints and variables are stored as: first A's, then B's. /// /// The product tableau has row layout: /// A's redundant rows, B's redundant rows, A's other rows, B's other rows. /// /// It has column layout: /// denominator, constant, A's columns, B's columns. Simplex Simplex::makeProduct(const Simplex &a, const Simplex &b) { … } std::optional<SmallVector<Fraction, 8>> Simplex::getRationalSample() const { … } void LexSimplexBase::addInequality(ArrayRef<DynamicAPInt> coeffs) { … } MaybeOptimum<SmallVector<Fraction, 8>> LexSimplex::getRationalSample() const { … } std::optional<SmallVector<DynamicAPInt, 8>> Simplex::getSamplePointIfIntegral() const { … } /// Given a simplex for a polytope, construct a new simplex whose variables are /// identified with a pair of points (x, y) in the original polytope. Supports /// some operations needed for generalized basis reduction. In what follows, /// dotProduct(x, y) = x_1 * y_1 + x_2 * y_2 + ... x_n * y_n where n is the /// dimension of the original polytope. /// /// This supports adding equality constraints dotProduct(dir, x - y) == 0. It /// also supports rolling back this addition, by maintaining a snapshot stack /// that contains a snapshot of the Simplex's state for each equality, just /// before that equality was added. class presburger::GBRSimplex { … }; /// Reduce the basis to try and find a direction in which the polytope is /// "thin". This only works for bounded polytopes. /// /// This is an implementation of the algorithm described in the paper /// "An Implementation of Generalized Basis Reduction for Integer Programming" /// by W. Cook, T. Rutherford, H. E. Scarf, D. Shallcross. /// /// Let b_{level}, b_{level + 1}, ... b_n be the current basis. /// Let width_i(v) = max <v, x - y> where x and y are points in the original /// polytope such that <b_j, x - y> = 0 is satisfied for all level <= j < i. /// /// In every iteration, we first replace b_{i+1} with b_{i+1} + u*b_i, where u /// is the integer such that width_i(b_{i+1} + u*b_i) is minimized. Let dual_i /// be the dual variable associated with the constraint <b_i, x - y> = 0 when /// computing width_{i+1}(b_{i+1}). It can be shown that dual_i is the /// minimizing value of u, if it were allowed to be fractional. Due to /// convexity, the minimizing integer value is either floor(dual_i) or /// ceil(dual_i), so we just need to check which of these gives a lower /// width_{i+1} value. If dual_i turned out to be an integer, then u = dual_i. /// /// Now if width_i(b_{i+1}) < 0.75 * width_i(b_i), we swap b_i and (the new) /// b_{i + 1} and decrement i (unless i = level, in which case we stay at the /// same i). Otherwise, we increment i. /// /// We keep f values and duals cached and invalidate them when necessary. /// Whenever possible, we use them instead of recomputing them. We implement the /// algorithm as follows. /// /// In an iteration at i we need to compute: /// a) width_i(b_{i + 1}) /// b) width_i(b_i) /// c) the integer u that minimizes width_i(b_{i + 1} + u*b_i) /// /// If width_i(b_i) is not already cached, we compute it. /// /// If the duals are not already cached, we compute width_{i+1}(b_{i+1}) and /// store the duals from this computation. /// /// We call updateBasisWithUAndGetFCandidate, which finds the minimizing value /// of u as explained before, caches the duals from this computation, sets /// b_{i+1} to b_{i+1} + u*b_i, and returns the new value of width_i(b_{i+1}). /// /// Now if width_i(b_{i+1}) < 0.75 * width_i(b_i), we swap b_i and b_{i+1} and /// decrement i, resulting in the basis /// ... b_{i - 1}, b_{i + 1} + u*b_i, b_i, b_{i+2}, ... /// with corresponding f values /// ... width_{i-1}(b_{i-1}), width_i(b_{i+1} + u*b_i), width_{i+1}(b_i), ... /// The values up to i - 1 remain unchanged. We have just gotten the middle /// value from updateBasisWithUAndGetFCandidate, so we can update that in the /// cache. The value at width_{i+1}(b_i) is unknown, so we evict this value from /// the cache. The iteration after decrementing needs exactly the duals from the /// computation of width_i(b_{i + 1} + u*b_i), so we keep these in the cache. /// /// When incrementing i, no cached f values get invalidated. However, the cached /// duals do get invalidated as the duals for the higher levels are different. void Simplex::reduceBasis(IntMatrix &basis, unsigned level) { … } /// Search for an integer sample point using a branch and bound algorithm. /// /// Each row in the basis matrix is a vector, and the set of basis vectors /// should span the space. Initially this is the identity matrix, /// i.e., the basis vectors are just the variables. /// /// In every level, a value is assigned to the level-th basis vector, as /// follows. Compute the minimum and maximum rational values of this direction. /// If only one integer point lies in this range, constrain the variable to /// have this value and recurse to the next variable. /// /// If the range has multiple values, perform generalized basis reduction via /// reduceBasis and then compute the bounds again. Now we try constraining /// this direction in the first value in this range and "recurse" to the next /// level. If we fail to find a sample, we try assigning the direction the next /// value in this range, and so on. /// /// If no integer sample is found from any of the assignments, or if the range /// contains no integer value, then of course the polytope is empty for the /// current assignment of the values in previous levels, so we return to /// the previous level. /// /// If we reach the last level where all the variables have been assigned values /// already, then we simply return the current sample point if it is integral, /// and go back to the previous level otherwise. /// /// To avoid potentially arbitrarily large recursion depths leading to stack /// overflows, this algorithm is implemented iteratively. std::optional<SmallVector<DynamicAPInt, 8>> Simplex::findIntegerSample() { … } /// Compute the minimum and maximum integer values the expression can take. We /// compute each separately. std::pair<MaybeOptimum<DynamicAPInt>, MaybeOptimum<DynamicAPInt>> Simplex::computeIntegerBounds(ArrayRef<DynamicAPInt> coeffs) { … } bool Simplex::isFlatAlong(ArrayRef<DynamicAPInt> coeffs) { … } void SimplexBase::print(raw_ostream &os) const { … } void SimplexBase::dump() const { … } bool Simplex::isRationalSubsetOf(const IntegerRelation &rel) { … } /// Returns the type of the inequality with coefficients `coeffs`. /// Possible types are: /// Redundant The inequality is satisfied by all points in the polytope /// Cut The inequality is satisfied by some points, but not by others /// Separate The inequality is not satisfied by any point /// /// Internally, this computes the minimum and the maximum the inequality with /// coefficients `coeffs` can take. If the minimum is >= 0, the inequality holds /// for all points in the polytope, so it is redundant. If the minimum is <= 0 /// and the maximum is >= 0, the points in between the minimum and the /// inequality do not satisfy it, the points in between the inequality and the /// maximum satisfy it. Hence, it is a cut inequality. If both are < 0, no /// points of the polytope satisfy the inequality, which means it is a separate /// inequality. Simplex::IneqType Simplex::findIneqType(ArrayRef<DynamicAPInt> coeffs) { … } /// Checks whether the type of the inequality with coefficients `coeffs` /// is Redundant. bool Simplex::isRedundantInequality(ArrayRef<DynamicAPInt> coeffs) { … } /// Check whether the equality given by `coeffs == 0` is redundant given /// the existing constraints. This is redundant when `coeffs` is already /// always zero under the existing constraints. `coeffs` is always zero /// when the minimum and maximum value that `coeffs` can take are both zero. bool Simplex::isRedundantEquality(ArrayRef<DynamicAPInt> coeffs) { … }
Codebase Browser
llvm