diff --git a/mlir/include/mlir/Analysis/AffineStructures.h b/mlir/include/mlir/Analysis/AffineStructures.h --- a/mlir/include/mlir/Analysis/AffineStructures.h +++ b/mlir/include/mlir/Analysis/AffineStructures.h @@ -13,6 +13,7 @@ #ifndef MLIR_ANALYSIS_AFFINE_STRUCTURES_H #define MLIR_ANALYSIS_AFFINE_STRUCTURES_H +#include "mlir/Analysis/Presburger/Matrix.h" #include "mlir/IR/AffineExpr.h" #include "mlir/IR/OpDefinition.h" #include "mlir/Support/LogicalResult.h" @@ -153,6 +154,12 @@ /// false if a solution exists or all tests were inconclusive. bool isIntegerEmpty() const; + // Returns a matrix where each row is a vector along which the polytope is + // bounded. The span of the returned vectors is guaranteed to contain all + // such vectors. The returned vectors are NOT guaranteed to be linearly + // independent. This function should not be called on empty sets. + Matrix getBoundedDirections() const; + /// Find a sample point satisfying the constraints. This uses a branch and /// bound algorithm with generalized basis reduction, which always works if /// the set is bounded. This should not be called for unbounded sets. diff --git a/mlir/include/mlir/Analysis/LinearTransform.h b/mlir/include/mlir/Analysis/LinearTransform.h new file mode 100644 --- /dev/null +++ b/mlir/include/mlir/Analysis/LinearTransform.h @@ -0,0 +1,47 @@ +//===- LinearTransform.h - MLIR LinearTransform Class -----------*- C++ -*-===// +// +// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions. +// See https://llvm.org/LICENSE.txt for license information. +// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception +// +//===----------------------------------------------------------------------===// +// +// Support for linear transforms and applying them to FlatAffineConstraints. +// +//===----------------------------------------------------------------------===// + +#ifndef MLIR_ANALYSIS_LINEARTRANSFORM_H +#define MLIR_ANALYSIS_LINEARTRANSFORM_H + +#include "mlir/Analysis/AffineStructures.h" +#include "mlir/Analysis/Presburger/Matrix.h" +#include "llvm/ADT/SmallVector.h" + +namespace mlir { + +class LinearTransform { +public: + // Returns a linear transform T such that MT is M in column echelon form. + // Also returns the number of non-zero columns in MT. + // + // Specifically, T is such that in every column the first non-zero row is + // strictly below that of the previous column, and all columns which have only + // zeros are at the end. + static std::pair + makeTransformToColumnEchelon(Matrix m); + + // Post-multiply every constraint vector v in fac with this transform, say T, + // turning it into vT. + FlatAffineConstraints applyTo(const FlatAffineConstraints &fac); + +private: + // Post-multiply the given vector v with this transform, say T, returning vT. + SmallVector postMultiplyVector(ArrayRef v); + explicit LinearTransform(Matrix &&oMatrix); + explicit LinearTransform(const Matrix &oMatrix); + + Matrix matrix; +}; + +} // namespace mlir +#endif // MLIR_ANALYSIS_LINEARTRANSFORM_H diff --git a/mlir/include/mlir/Analysis/Presburger/Fraction.h b/mlir/include/mlir/Analysis/Presburger/Fraction.h --- a/mlir/include/mlir/Analysis/Presburger/Fraction.h +++ b/mlir/include/mlir/Analysis/Presburger/Fraction.h @@ -64,6 +64,8 @@ inline bool operator==(Fraction x, Fraction y) { return compare(x, y) == 0; } +inline bool operator!=(Fraction x, Fraction y) { return compare(x, y) != 0; } + inline bool operator>(Fraction x, Fraction y) { return compare(x, y) > 0; } inline bool operator>=(Fraction x, Fraction y) { return compare(x, y) >= 0; } diff --git a/mlir/include/mlir/Analysis/Presburger/Matrix.h b/mlir/include/mlir/Analysis/Presburger/Matrix.h --- a/mlir/include/mlir/Analysis/Presburger/Matrix.h +++ b/mlir/include/mlir/Analysis/Presburger/Matrix.h @@ -58,6 +58,12 @@ /// Add `scale` multiples of the source row to the target row. void addToRow(unsigned sourceRow, unsigned targetRow, int64_t scale); + /// Add `scale` multiples of the source column to the target column. + void addToColumn(unsigned sourceColumn, unsigned targetColumn, int64_t scale); + + /// Negate the specified column. + void negateColumn(unsigned column); + /// Resize the matrix to the specified dimensions. If a dimension is smaller, /// the values are truncated; if it is bigger, the new values are default /// initialized. diff --git a/mlir/include/mlir/Analysis/Presburger/Simplex.h b/mlir/include/mlir/Analysis/Presburger/Simplex.h --- a/mlir/include/mlir/Analysis/Presburger/Simplex.h +++ b/mlir/include/mlir/Analysis/Presburger/Simplex.h @@ -17,10 +17,12 @@ #include "mlir/Analysis/AffineStructures.h" #include "mlir/Analysis/Presburger/Fraction.h" #include "mlir/Analysis/Presburger/Matrix.h" +#include "mlir/IR/Location.h" #include "mlir/Support/LogicalResult.h" #include "llvm/ADT/ArrayRef.h" #include "llvm/ADT/Optional.h" #include "llvm/ADT/SmallVector.h" +#include "llvm/Support/StringSaver.h" #include "llvm/Support/raw_ostream.h" namespace mlir { @@ -84,7 +86,7 @@ /// /// The unknowns in row position are represented in terms of the basis unknowns. /// If the basis unknowns are u_1, u_2, ... u_m, and a row in the tableau is -/// d, c, a_1, a_2, ... a_m, this representats the unknown for that row as +/// d, c, a_1, a_2, ... a_m, this represents the unknown for that row as /// (c + a_1*u_1 + a_2*u_2 + ... + a_m*u_m)/d. In our running example, if the /// basis is the initial basis of x, y, then the constraint 1 + 2x + 3y >= 0 /// would be represented by the row [1, 1, 2, 3]. @@ -173,20 +175,25 @@ void intersectFlatAffineConstraints(const FlatAffineConstraints &fac); /// Compute the maximum or minimum value of the given row, depending on - /// direction. The specified row is never pivoted. + /// direction. The specified row is never pivoted. On return, the row may + /// have a negative sample value if the direction is down. /// - /// Returns a (num, den) pair denoting the optimum, or None if no - /// optimum exists, i.e., if the expression is unbounded in this direction. + /// Returns a Fraction denoting the optimum, or a null value if no optimum + /// exists, i.e., if the expression is unbounded in this direction. Optional computeRowOptimum(Direction direction, unsigned row); /// Compute the maximum or minimum value of the given expression, depending on - /// direction. + /// direction. Should not be called when the Simplex is empty. /// - /// Returns a (num, den) pair denoting the optimum, or a null value if no - /// optimum exists, i.e., if the expression is unbounded in this direction. + /// Returns a Fraction denoting the optimum, or a null value if no optimum + /// exists, i.e., if the expression is unbounded in this direction. Optional computeOptimum(Direction direction, ArrayRef coeffs); + /// Returns whether the perpendicular of the specified constraint is a + /// is a direction along which the polytope is bounded. + bool isBoundedAlongConstraint(unsigned constraintIndex); + /// Returns whether the specified constraint has been marked as redundant. /// Constraints are numbered from 0 starting at the first added inequality. /// Equalities are added as a pair of inequalities and so correspond to two @@ -299,6 +306,15 @@ /// sample value, false otherwise. LogicalResult restoreRow(Unknown &u); + /// Compute the maximum or minimum of the specified Unknown, depending on + /// direction. The specified unknown may be pivoted. If the unknown is + /// restricted, it will have a non-negative sample value on return. + /// Should not be called if the Simplex is empty. + /// + /// Returns a Fraction denoting the optimum, or a null value if no optimum + /// exists, i.e., if the expression is unbounded in this direction. + Optional computeOptimum(Direction direction, Unknown &u); + /// Mark the specified unknown redundant. This operation is added to the undo /// log and will be undone by rollbacks. The specified unknown must be in row /// orientation. diff --git a/mlir/lib/Analysis/AffineStructures.cpp b/mlir/lib/Analysis/AffineStructures.cpp --- a/mlir/lib/Analysis/AffineStructures.cpp +++ b/mlir/lib/Analysis/AffineStructures.cpp @@ -11,6 +11,7 @@ //===----------------------------------------------------------------------===// #include "mlir/Analysis/AffineStructures.h" +#include "mlir/Analysis/LinearTransform.h" #include "mlir/Analysis/Presburger/Simplex.h" #include "mlir/Dialect/Affine/IR/AffineOps.h" #include "mlir/Dialect/Affine/IR/AffineValueMap.h" @@ -20,6 +21,7 @@ #include "mlir/Support/LLVM.h" #include "mlir/Support/MathExtras.h" #include "llvm/ADT/SmallPtrSet.h" +#include "llvm/ADT/SmallVector.h" #include "llvm/Support/Debug.h" #include "llvm/Support/raw_ostream.h" @@ -1034,21 +1036,152 @@ return false; } -// First, try the GCD test heuristic. +// Returns a matrix where each row is a vector along which the polytope is +// bounded. The span of the returned vectors is guaranteed to contain all +// such vectors. The returned vectors are NOT guaranteed to be linearly +// independent. This function should not be called on empty sets. // -// If that doesn't find the set empty, check if the set is unbounded. If it is, -// we cannot use the GBR algorithm and we conservatively return false. -// -// If the set is bounded, we use the complete emptiness check for this case -// provided by Simplex::findIntegerSample(), which gives a definitive answer. +// It is sufficient to check the perpendiculars of the constraints, as the set +// of perpendiculars which are bounded must span all bounded directions. +Matrix FlatAffineConstraints::getBoundedDirections() const { + // Note that it is necessary to add the equalities too (which the constructor + // does) even though we don't need to check if they are bounded; whether an + // inequality is bounded or not depends on what other constraints, including + // equalities, are present. + Simplex simplex(*this); + + assert(!simplex.isEmpty() && "It is not meaningful to ask whether a " + "direction is bounded in an empty set."); + + SmallVector boundedIneqs; + // The constructor adds the inequalities to the simplex first, so this + // processes all the inequalities. + for (unsigned i = 0, e = getNumInequalities(); i < e; ++i) { + if (simplex.isBoundedAlongConstraint(i)) + boundedIneqs.push_back(i); + } + + // The direction vector is given by the coefficients and does not include the + // constant term, so the matrix has one fewer column. + const unsigned dirsNumCols = getNumCols() - 1; + Matrix dirs(boundedIneqs.size() + getNumEqualities(), dirsNumCols); + + // Copy the bounded inequalities. + unsigned row = 0; + for (unsigned i : boundedIneqs) { + for (unsigned col = 0; col < dirsNumCols; ++col) + dirs(row, col) = atIneq(i, col); + ++row; + } + + // Copy the equalities. All the equalities' perpendiculars are bounded. + for (unsigned i = 0, e = getNumEqualities(); i < e; ++i) { + for (unsigned col = 0; col < dirsNumCols; ++col) + dirs(row, col) = atEq(i, col); + ++row; + } + + return dirs; +} + +bool eqInvolvesSuffixDims(const FlatAffineConstraints &fac, unsigned eqIndex, + unsigned numDims) { + for (unsigned e = fac.getNumDimIds(), j = e - numDims; j < e; ++j) + if (fac.atEq(eqIndex, j) != 0) + return true; + return false; +} +bool ineqInvolvesSuffixDims(const FlatAffineConstraints &fac, + unsigned ineqIndex, unsigned numDims) { + for (unsigned e = fac.getNumDimIds(), j = e - numDims; j < e; ++j) + if (fac.atIneq(ineqIndex, j) != 0) + return true; + return false; +} + +void removeConstraintsInvolvingSuffixDims(FlatAffineConstraints &fac, + unsigned unboundedDims) { + // We iterate backwards so that whether we remove constraint i - 1 or not, the + // next constraint to be tested is always i - 2. + for (unsigned i = fac.getNumEqualities(); i > 0; i--) + if (eqInvolvesSuffixDims(fac, i - 1, unboundedDims)) + fac.removeEquality(i - 1); + for (unsigned i = fac.getNumInequalities(); i > 0; i--) + if (ineqInvolvesSuffixDims(fac, i - 1, unboundedDims)) + fac.removeInequality(i - 1); +} + +/// Let this set be S. If S is bounded then we directly call into the GBR +/// sampling algorithm. Otherwise, there are some unbounded directions, i.e., +/// vectors v such that S extends to infininty along v or -v. In this case we +/// use an algorithm described in the integer set library (isl) manual and used +/// by the isl_set_sample function in that library. The algorithm is: +/// +/// 1) Apply a unimodular transform T to S to obtain S*T, such that all +/// dimensions in which S*T is bounded lie in the linear span of a prefix of the +/// dimensions. +/// +/// 2) Construct a set transformedSet by removing all constraints that involve +/// the unbounded dimensions and also deleting the unbounded dimensions. Note +/// that this is a bounded set. +/// +/// 3) Check if transformedSet is empty using the GBR sampling algorithm. +/// +/// 4) return S is empty iff transformedSet is empty. +/// +/// Since T is unimodular, a vector v is a solution to S*T iff T*v is a +/// solution to S. The following is a sketch of a proof that S*T is empty +/// iff transformedSet is empty: +/// +/// If transformedSet is empty, then S*T is certainly empty since transformedSet +/// was obtained by removing constraints and deleting dimensions from S*T. +/// +/// If transformedSet contains a sample, consider the set C obtained by +/// substituting the sample for the bounded dimensions of S*T. All the +/// constraints of S*T that did not involve unbounded dimensions are +/// satisfied by this substitution. +/// +/// In step 1, all dimensions in the linear span of the dimensions outside the +/// prefix are unbounded in S*T. Substituting values for the bounded dimensions +/// cannot makes these dimensions bounded, and these are the only remaining +/// dimensions in C, so C is unbounded along every vector. C is hence a +/// full-dimensional cone and therefore always contains an integer point, which +/// we can then substitute to get a full solution to S*T. bool FlatAffineConstraints::isIntegerEmpty() const { + // First, try the GCD test heuristic. if (isEmptyByGCDTest()) return true; Simplex simplex(*this); - if (simplex.isUnbounded()) - return false; - return !simplex.findIntegerSample().hasValue(); + if (simplex.isEmpty()) + return true; + + // For a bounded set, we directly call into the GBR sampling algorithm. + if (!simplex.isUnbounded()) + return !simplex.findIntegerSample().hasValue(); + + // The set is unbounded. We cannot directly use the GBR algorithm. + // + // m is a matrix containing, in each row, a vector in which S is + // bounded, such that the linear span of all these dimensions contains all + // bounded dimensions in S. + Matrix m = getBoundedDirections(); + // In column echelon form, each row of m occupies only the first rank(m) + // columns and has zeros on the other columns. The transform T that brings S + // to column echelon form is unimodular as well, so this is a suitable + // transform to use in step 1 of the algorithm. + std::pair result = + LinearTransform::makeTransformToColumnEchelon(std::move(m)); + FlatAffineConstraints transformedSet = result.second.applyTo(*this); + + unsigned numBoundedDims = result.first; + unsigned numUnboundedDims = getNumIds() - numBoundedDims; + removeConstraintsInvolvingSuffixDims(transformedSet, numUnboundedDims); + + // Remove all the unbounded dimensions. + transformedSet.removeIdRange(numBoundedDims, transformedSet.getNumIds()); + + return !Simplex(transformedSet).findIntegerSample().hasValue(); } Optional> diff --git a/mlir/lib/Analysis/CMakeLists.txt b/mlir/lib/Analysis/CMakeLists.txt --- a/mlir/lib/Analysis/CMakeLists.txt +++ b/mlir/lib/Analysis/CMakeLists.txt @@ -3,6 +3,7 @@ AffineStructures.cpp BufferAliasAnalysis.cpp CallGraph.cpp + LinearTransform.cpp Liveness.cpp LoopAnalysis.cpp NestedMatcher.cpp @@ -36,6 +37,7 @@ add_mlir_library(MLIRLoopAnalysis AffineAnalysis.cpp AffineStructures.cpp + LinearTransform.cpp LoopAnalysis.cpp NestedMatcher.cpp PresburgerSet.cpp diff --git a/mlir/lib/Analysis/LinearTransform.cpp b/mlir/lib/Analysis/LinearTransform.cpp new file mode 100644 --- /dev/null +++ b/mlir/lib/Analysis/LinearTransform.cpp @@ -0,0 +1,157 @@ +//===- LinearTransform.cpp - MLIR LinearTransform Class -------------------===// +// +// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions. +// See https://llvm.org/LICENSE.txt for license information. +// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception +// +//===----------------------------------------------------------------------===// + +#include "mlir/Analysis/LinearTransform.h" +#include "mlir/Analysis/AffineStructures.h" + +namespace mlir { + +LinearTransform::LinearTransform(Matrix &&oMatrix) : matrix(oMatrix) {} +LinearTransform::LinearTransform(const Matrix &oMatrix) : matrix(oMatrix) {} + +// Set M(row, targetCol) to its remainder on division by M(row, sourceCol) +// by subtracting from column targetCol an appropriate integer multiple of +// sourceCol. This brings M(row, targetCol) to the range [0, M(row, sourceCol)). +// Apply the same column operation to otherMatrix, with the same integer +// multiple. +static void modEntryColumnOperation(Matrix &m, unsigned row, unsigned sourceCol, + unsigned targetCol, Matrix &otherMatrix) { + assert(m(row, sourceCol) != 0 && "Cannot divide by zero!"); + assert((m(row, sourceCol) > 0 && m(row, targetCol) > 0) && + "Operands must be positive!"); + int64_t ratio = m(row, targetCol) / m(row, sourceCol); + m.addToColumn(sourceCol, targetCol, -ratio); + otherMatrix.addToColumn(sourceCol, targetCol, -ratio); +} + +std::pair +LinearTransform::makeTransformToColumnEchelon(Matrix m) { + // We start with an identity result matrix and perform operations on m + // until m is in column echelon form. We apply the same sequence of operations + // on resultMatrix to obtain a transform that takes m to column echelon + // form. + Matrix resultMatrix = Matrix::identity(m.getNumColumns()); + + unsigned echelonCol = 0; + // Invariant: in all rows above row, all columns from echelonCol onwards + // are all zero elements. In an iteration, if the curent row has any non-zero + // elements echelonCol onwards, we bring one to echelonCol and use it to + // make all elements echelonCol + 1 onwards zero. + for (unsigned row = 0; row < m.getNumRows(); ++row) { + // Search row for a non-empty entry, starting at echelonCol. + unsigned nonZeroCol = echelonCol; + for (unsigned e = m.getNumColumns(); nonZeroCol < e; ++nonZeroCol) { + if (m(row, nonZeroCol) == 0) + continue; + break; + } + + // Continue to the next row with the same echelonCol if this row is all + // zeros from echelonCol onwards. + if (nonZeroCol == m.getNumColumns()) + continue; + + // Bring the non-zero column to echelonCol. This doesn't affect rows + // above since they are all zero at these columns. + if (nonZeroCol != echelonCol) { + m.swapColumns(nonZeroCol, echelonCol); + resultMatrix.swapColumns(nonZeroCol, echelonCol); + } + + // Make m(row, echelonCol) non-negative. + if (m(row, echelonCol) < 0) { + m.negateColumn(echelonCol); + resultMatrix.negateColumn(echelonCol); + } + + // Make all the entries in row after echelonCol zero. + for (unsigned i = echelonCol + 1, e = m.getNumColumns(); i < e; ++i) { + // We make m(row, i) non-negative, and then apply the Euclidean GCD + // algorithm to (row, i) and (row, echelonCol). At the end, one of them + // has value equal to the gcd of the two entries, and the other is zero. + + if (m(row, i) < 0) { + m.negateColumn(i); + resultMatrix.negateColumn(i); + } + + unsigned targetCol = i, sourceCol = echelonCol; + // At every step, we set m(row, targetCol) %= m(row, sourceCol), and + // swap the indices sourceCol and targetCol. (not the columns themselves) + // This modulo is implemented as a subtraction + // m(row, targetCol) -= quotient * m(row, sourceCol), + // where quotient = floor(m(row, targetCol) / m(row, sourceCol)), + // which brings m(row, targetCol) to the range [0, m(row, sourceCol)). + // + // We are only allowed column operations; we perform the above + // for every row, i.e., the above subtraction is done as a column + // operation. This does not affect any rows above us since they are + // guaranteed to be zero at these columns. + while (m(row, targetCol) != 0 && m(row, sourceCol) != 0) { + modEntryColumnOperation(m, row, sourceCol, targetCol, resultMatrix); + std::swap(targetCol, sourceCol); + } + + // One of (row, echelonCol) and (row, i) is zero and the other is the gcd. + // Make it so that (row, echelonCol) holds the non-zero value. + if (m(row, echelonCol) == 0) { + m.swapColumns(i, echelonCol); + resultMatrix.swapColumns(i, echelonCol); + } + } + + ++echelonCol; + } + + return {echelonCol, LinearTransform(std::move(resultMatrix))}; +} + +SmallVector +LinearTransform::postMultiplyVector(ArrayRef v) { + assert(v.size() == matrix.getNumRows() && + "vector dimension should be matrix output dimension"); + + SmallVector result; + result.reserve(v.size()); + for (unsigned col = 0, e = matrix.getNumColumns(); col < e; ++col) { + int64_t elem = 0; + for (unsigned i = 0, e = matrix.getNumRows(); i < e; ++i) + elem += v[i] * matrix(i, col); + result.push_back(elem); + } + return result; +} + +FlatAffineConstraints +LinearTransform::applyTo(const FlatAffineConstraints &fac) { + FlatAffineConstraints result(fac.getNumDimIds()); + + for (unsigned i = 0, e = fac.getNumEqualities(); i < e; ++i) { + ArrayRef eq = fac.getEquality(i); + + int64_t c = eq.back(); + + SmallVector newEq = postMultiplyVector(eq.drop_back()); + newEq.push_back(c); + result.addEquality(newEq); + } + + for (unsigned i = 0, e = fac.getNumInequalities(); i < e; ++i) { + ArrayRef ineq = fac.getInequality(i); + + int64_t c = ineq.back(); + + SmallVector newIneq = postMultiplyVector(ineq.drop_back()); + newIneq.push_back(c); + result.addInequality(newIneq); + } + + return result; +} + +} // namespace mlir \ No newline at end of file diff --git a/mlir/lib/Analysis/Presburger/CMakeLists.txt b/mlir/lib/Analysis/Presburger/CMakeLists.txt --- a/mlir/lib/Analysis/Presburger/CMakeLists.txt +++ b/mlir/lib/Analysis/Presburger/CMakeLists.txt @@ -1,4 +1,4 @@ add_mlir_library(MLIRPresburger Simplex.cpp Matrix.cpp - ) \ No newline at end of file + ) diff --git a/mlir/lib/Analysis/Presburger/Matrix.cpp b/mlir/lib/Analysis/Presburger/Matrix.cpp --- a/mlir/lib/Analysis/Presburger/Matrix.cpp +++ b/mlir/lib/Analysis/Presburger/Matrix.cpp @@ -79,6 +79,19 @@ return; } +void Matrix::addToColumn(unsigned sourceColumn, unsigned targetColumn, + int64_t scale) { + if (scale == 0) + return; + for (unsigned row = 0, e = getNumRows(); row < e; ++row) + at(row, targetColumn) += scale * at(row, sourceColumn); +} + +void Matrix::negateColumn(unsigned column) { + for (unsigned row = 0, e = getNumRows(); row < e; ++row) + at(row, column) = -at(row, column); +} + void Matrix::print(raw_ostream &os) const { for (unsigned row = 0; row < nRows; ++row) { for (unsigned column = 0; column < nColumns; ++column) diff --git a/mlir/lib/Analysis/Presburger/Simplex.cpp b/mlir/lib/Analysis/Presburger/Simplex.cpp --- a/mlir/lib/Analysis/Presburger/Simplex.cpp +++ b/mlir/lib/Analysis/Presburger/Simplex.cpp @@ -9,6 +9,7 @@ #include "mlir/Analysis/Presburger/Simplex.h" #include "mlir/Analysis/Presburger/Matrix.h" #include "mlir/Support/MathExtras.h" +#include "llvm/ADT/Optional.h" namespace mlir { using Direction = Simplex::Direction; @@ -482,7 +483,7 @@ /// or None if it is unbounded. Optional Simplex::computeOptimum(Direction direction, ArrayRef coeffs) { - assert(!empty && "Tableau should not be empty"); + assert(!empty && "Simplex should not be empty"); unsigned snapshot = getSnapshot(); unsigned conIndex = addRow(coeffs); @@ -492,6 +493,34 @@ return optimum; } +Optional Simplex::computeOptimum(Direction direction, Unknown &u) { + assert(!empty && "Simplex should not be empty!"); + if (u.orientation == Orientation::Column) { + unsigned column = u.pos; + Optional pivotRow = findPivotRow({}, direction, column); + // If no pivot is returned, the constraint is unbounded in the specified + // direction. + if (!pivotRow) + return {}; + pivot(*pivotRow, column); + } + + unsigned row = u.pos; + Optional optimum = computeRowOptimum(direction, row); + if (u.restricted && direction == Direction::Down && + (!optimum || *optimum < Fraction(0, 1))) + restoreRow(u); + return optimum; +} + +bool Simplex::isBoundedAlongConstraint(unsigned constraintIndex) { + assert(!empty && "It is not meaningful to ask whether a direction is bounded " + "in an empty set."); + // The constraint's perpendicular is already bounded below, since it is a + // constraint. If it is also bounded above, we can return true. + return computeOptimum(Direction::Up, con[constraintIndex]).hasValue(); +} + /// Redundant constraints are those that are in row orientation and lie in /// rows 0 to nRedundant - 1. bool Simplex::isMarkedRedundant(unsigned constraintIndex) const { diff --git a/mlir/unittests/Analysis/AffineStructuresTest.cpp b/mlir/unittests/Analysis/AffineStructuresTest.cpp --- a/mlir/unittests/Analysis/AffineStructuresTest.cpp +++ b/mlir/unittests/Analysis/AffineStructuresTest.cpp @@ -15,22 +15,36 @@ namespace mlir { -/// If 'hasValue' is true, check that findIntegerSample returns a valid sample +enum class TestFunction { Sample, Empty }; + +/// If fn is TestFunction::Sample (default): +/// If hasSample is true, check that findIntegerSample returns a valid sample /// for the FlatAffineConstraints fac. +/// If hasSample is false, check that findIntegerSample returns None. /// -/// If hasValue is false, check that findIntegerSample does not return None. -static void checkSample(bool hasValue, const FlatAffineConstraints &fac) { - Optional> maybeSample = fac.findIntegerSample(); - if (!hasValue) { - EXPECT_FALSE(maybeSample.hasValue()); - if (maybeSample.hasValue()) { - for (auto x : *maybeSample) - llvm::errs() << x << ' '; - llvm::errs() << '\n'; +/// If fn is TestFunction::Empty, check that isIntegerEmpty returns the +/// opposite of hasSample. +static void checkSample(bool hasSample, const FlatAffineConstraints &fac, + TestFunction fn = TestFunction::Sample) { + Optional> maybeSample; + switch (fn) { + case TestFunction::Sample: + maybeSample = fac.findIntegerSample(); + if (!hasSample) { + EXPECT_FALSE(maybeSample.hasValue()); + if (maybeSample.hasValue()) { + for (auto x : *maybeSample) + llvm::errs() << x << ' '; + llvm::errs() << '\n'; + } + } else { + ASSERT_TRUE(maybeSample.hasValue()); + EXPECT_TRUE(fac.containsPoint(*maybeSample)); } - } else { - ASSERT_TRUE(maybeSample.hasValue()); - EXPECT_TRUE(fac.containsPoint(*maybeSample)); + break; + case TestFunction::Empty: + EXPECT_EQ(!hasSample, fac.isIntegerEmpty()); + break; } } @@ -52,9 +66,11 @@ /// orderings may cause the algorithm to proceed differently. At least some of ///.these permutations should make it past the heuristics and test the /// implementation of the GBR algorithm itself. -static void checkPermutationsSample(bool hasValue, unsigned nDim, +/// Use TestFunction fn to test. +static void checkPermutationsSample(bool hasSample, unsigned nDim, ArrayRef> ineqs, - ArrayRef> eqs) { + ArrayRef> eqs, + TestFunction fn = TestFunction::Sample) { SmallVector perm(nDim); std::iota(perm.begin(), perm.end(), 0); auto permute = [&perm](ArrayRef coeffs) { @@ -71,8 +87,8 @@ for (const auto &eq : eqs) permutedEqs.push_back(permute(eq)); - checkSample(hasValue, - makeFACFromConstraints(nDim, permutedIneqs, permutedEqs)); + checkSample(hasSample, + makeFACFromConstraints(nDim, permutedIneqs, permutedEqs), fn); } while (std::next_permutation(perm.begin(), perm.end())); } @@ -206,19 +222,158 @@ EXPECT_FALSE( makeFACFromConstraints(1, {{5, -1}, {-5, 9}}, {}).isIntegerEmpty()); - // An unbounded set, which isIntegerEmpty should detect as unbounded and - // return without calling findIntegerSample. + // Unbounded sets. + EXPECT_TRUE(makeFACFromConstraints(3, + { + {2, 0, 0, -1}, // 2x >= 1 + {-2, 0, 0, 1}, // 2x <= 1 + {0, 2, 0, -1}, // 2y >= 1 + {0, -2, 0, 1}, // 2y <= 1 + {0, 0, 2, -1}, // 2z >= 1 + }, + {}) + .isIntegerEmpty()); + EXPECT_FALSE(makeFACFromConstraints(3, { - {2, 0, 0, -1}, - {-2, 0, 0, 1}, - {0, 2, 0, -1}, - {0, -2, 0, 1}, - {0, 0, 2, -1}, + {2, 0, 0, -1}, // 2x >= 1 + {-3, 0, 0, 3}, // 3x <= 3 + {0, 0, 5, -6}, // 5z >= 6 + {0, 0, -7, 17}, // 7z <= 17 + {0, 3, 0, -2}, // 3y >= 2 }, {}) .isIntegerEmpty()); + // 2D cone with apex at (10000, 10000) and + // edges passing through (1/3, 0) and (2/3, 0). + EXPECT_FALSE( + makeFACFromConstraints( + 2, {{300000, -299999, -100000}, {-300000, 299998, 200000}}, {}) + .isIntegerEmpty()); + + // Cartesian product of a tetrahedron and a 2D cone. + // The tetrahedron has vertices at + // (1/3, 0, 0), (2/3, 0, 0), (2/3, 0, 10000), and (10000, 10000, 10000). + // The first three points form a triangular base on the xz plane with the + // apex at the fourth point, which is the only integer point. + // The cone has apex at (10000, 10000) and + // edges passing through (1/3, 0) and (2/3, 0). + checkPermutationsSample( + true /* not empty */, 5, + { + // Tetrahedron contraints: + {0, 1, 0, 0, 0, 0}, // y >= 0 + {0, -1, 1, 0, 0, 0}, // z >= y + // -300000x + 299998y + 100000 + z <= 0. + {300000, -299998, -1, 0, 0, -100000}, + // -150000x + 149999y + 100000 >= 0. + {-150000, 149999, 0, 0, 0, 100000}, + + // Triangle constraints: + // 300000p - 299999q >= 100000 + {0, 0, 0, 300000, -299999, -100000}, + // -300000p + 299998q + 200000 >= 0 + {0, 0, 0, -300000, 299998, 200000}, + }, + {}, TestFunction::Empty); + + // Cartesian product of same tetrahedron as above and {(p, q) : 1/3 <= p <= + // 2/3}. Since the second set is empty, the whole set is too. + checkPermutationsSample( + false /* empty */, 5, + { + // Tetrahedron contraints: + {0, 1, 0, 0, 0, 0}, // y >= 0 + {0, -1, 1, 0, 0, 0}, // z >= y + // -300000x + 299998y + 100000 + z <= 0. + {300000, -299998, -1, 0, 0, -100000}, + // -150000x + 149999y + 100000 >= 0. + {-150000, 149999, 0, 0, 0, 100000}, + + // Second set constraints: + // 3p >= 1 + {0, 0, 0, 3, 0, -1}, + // 3p <= 2 + {0, 0, 0, -3, 0, 2}, + }, + {}, TestFunction::Empty); + + // Cartesian product of same tetrahedron as above and + // {(p, q, r) : 1 <= p <= 2 and p = 3q + 3r}. + // Since the second set is empty, the whole set is too. + checkPermutationsSample( + false /* empty */, 5, + { + // Tetrahedron contraints: + {0, 1, 0, 0, 0, 0, 0}, // y >= 0 + {0, -1, 1, 0, 0, 0, 0}, // z >= y + // -300000x + 299998y + 100000 + z <= 0. + {300000, -299998, -1, 0, 0, 0, -100000}, + // -150000x + 149999y + 100000 >= 0. + {-150000, 149999, 0, 0, 0, 0, 100000}, + + // Second set constraints: + // p >= 1 + {0, 0, 0, 1, 0, 0, -1}, + // p <= 2 + {0, 0, 0, -1, 0, 0, 2}, + }, + { + {0, 0, 0, 1, -3, -3, 0}, // p = 3q + 3r + }, + TestFunction::Empty); + + // Cartesian product of a tetrahedron and a 2D cone. + // The tetrahedron is empty and has vertices at + // (1/3, 0, 0), (2/3, 0, 0), (2/3, 0, 100), and (100, 100 - 1/3, 100). + // The cone has apex at (10000, 10000) and + // edges passing through (1/3, 0) and (2/3, 0). + // Since the tetrahedron is empty, the Cartesian product is too. + checkPermutationsSample(false /* empty */, 5, + { + // Tetrahedron contraints: + {0, 1, 0, 0, 0, 0}, + {0, -300, 299, 0, 0, 0}, + {300 * 299, -89400, -299, 0, 0, -100 * 299}, + {-897, 894, 0, 0, 0, 598}, + + // Triangle constraints: + // 300000p - 299999q >= 100000 + {0, 0, 0, 300000, -299999, -100000}, + // -300000p + 299998q + 200000 >= 0 + {0, 0, 0, -300000, 299998, 200000}, + }, + {}, TestFunction::Empty); + + // Cartesian product of same tetrahedron as above and + // {(p, q) : 1/3 <= p <= 2/3}. + checkPermutationsSample(false /* empty */, 5, + { + // Tetrahedron contraints: + {0, 1, 0, 0, 0, 0}, + {0, -300, 299, 0, 0, 0}, + {300 * 299, -89400, -299, 0, 0, -100 * 299}, + {-897, 894, 0, 0, 0, 598}, + + // Second set constraints: + // 3p >= 1 + {0, 0, 0, 3, 0, -1}, + // 3p <= 2 + {0, 0, 0, -3, 0, 2}, + }, + {}, TestFunction::Empty); + + EXPECT_FALSE(makeFACFromConstraints(3, + { + {2, 0, 0, -1}, // 2x >= 1 + }, + {{ + {1, -1, 0, -1}, // y = x - 1 + {0, 1, -1, 0}, // z = y + }}) + .isIntegerEmpty()); + // FlatAffineConstraints::isEmpty() does not detect the following sets to be // empty.