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