Support emptiness checks for unbounded FlatAffineConstraints.

[arjunp]

Summary

With this, we have complete support for emptiness checks. This also paves the way for future support to check if two FlatAffineConstraints are equal.

[arjunp]

Added attribution for the algorithm.

[ftynse]

It looks like something went wrong when updating the second revision and the diff now only contains the comment change. Note that if you use arc diff on a branch of multiple commits, you must specify the commit that the branch is based on rather than HEAD^, which is merely the previous-to-head commit.

[arjunp]

Trying to fix the patch.

[arjunp]

Seems to be fixed now.

[arjunp]

Resolved clang-tidy issue.

[ftynse]

mlir/include/mlir/Analysis/Presburger/LinearTransform.h

31

Naming nit: a.apply(b) reads like applying b to a, but this function does the inverse. applyTo?

mlir/lib/Analysis/AffineStructures.cpp

1079

Shouldn't this be atEq?

1103

Let's not commit commented-out code

1106–1107

Or you can do

for (unsigned i = fac.getNumEqualities(); i >0; --i) 
  if (eqInvolvesSuffixDims(fac, i - 1, unboundedDims))
    fac.removeEquality(i - 1);

that is, replace uses of i with i - 1 (for long loops, just define a variable), and avoid implicit sign conversion.

1185

You may want to call removeIdRange instead. This will avoid the loop together with its signedness magic and a lot of copies that this does because of removing numBoundedDims instead of i.

mlir/lib/Analysis/Presburger/LinearTransform.cpp

14

I suppose you wanted Matrix &&oMatrix instead. Otherwise, this makes a copy of a matrix and than moves out of that copy. If you do need a copy, consider const Matrix & instead.

16

I don't know how to understand "remainder of division by [0, M(row, sourceCol))". Did you rather mean [0, M(row, sourceCol) be the range of values the remainder can take?

23

Please only use auto when it improves readability, e.g. for long iterator types or where type is obvious from context such as cast.

23–24

Note that machine integer division is rounding to zero, as opposed to Euclidean division that rounds to negative infinity. So -4 / 3 gives -1 with machine division and -2 with Euclidean division. Usually, in affine world, I'd expect the latter. If it is used further to compute the remainder, machine division will yield -4 - (-1 * 3) = -1, but one would normally expect the remainder to be non-negative -4 - (-2 * 3) = 2. Please verify that this code does what you expect it to do and, if so, add an appropriate comment (-1 = 2 mod 3 after all).

75

Except that the code in modEntryColumnOperation doesn't do floor().

110

Please reserve the space in the vector before pushing back in a loop.

mlir/unittests/Analysis/AffineStructuresTest.cpp

254

Please add a test that involves equalities, and a test that exercises Euclidean division when computing the row echelon form. Since these are unit tests, the latter can be easily done by testing the row echelon form computation directly.

[arjunp]

Address comments.

[arjunp]

Thanks for the review!

mlir/lib/Analysis/AffineStructures.cpp

1079

Sorry, yes. Fixed and added a test case.

mlir/lib/Analysis/Presburger/LinearTransform.cpp

16

Yes, fixed.

23–24

Yeah, my intention was to use Euclidean division.

We can't test for this because the Euclidean GCD algorithm is still correct and terminates with machine division (about efficiency, I am not sure). It is still correct with machine division since adding _any_ multiple of one operand to the other preserves their GCD. Termination:

Suppose we run the algorithm with integers x and y. If one of x or y is divides the other then the remainder issue doesn't come into play.

Otherwise, after one iteration we have -|x| < y < |x|.

If |y| > |x|, this means |y| decreases after this iteration. Otherwise, consider the next iteration after x, y are swapped, in which the operand of greater magnitude will be set to the remainder. That operand will decrease in magnitude then. So every two iterations one of the operands decreases in magnitude as long as neither is divisible by the other.

Once one operand becomes a multiple of the other, that one gets set to zero and we are done. To test that we use Euclidean division, we would have to directly test this static free function modEntryColumnOperation.

I feel it's easier to reason about the code if both operands are made positive, so I've changed it accordingly.

mlir/unittests/Analysis/AffineStructuresTest.cpp

254

Added tests involving equalities. See above regarding the unit test for LinearTransform. The column echelon form computation works with both Euclidean and machine division.

[andydavis1]

mlir/lib/Analysis/AffineStructures.cpp

1065

Might be worth capturing this fact in local variable: unsigned dirsNumCols = getNumCols() - 1; since it is use 3 times below.

[arjunp]

Addressed Andy's comment.

[bollu]

mlir/lib/Analysis/AffineStructures.cpp

1066

consider

const unsigned dirsNumCols = getNumCols() - 1;

to indicate that the value isn't changed throughout the computation.

mlir/lib/Analysis/Presburger/Matrix.cpp

88

nit: eliminate unnecessary return?

mlir/lib/Analysis/Presburger/Simplex.cpp

523

Consider:

return computeOptimum(Direction::Up, con[constraintIndex]).hasValue();

instead of:

if (Optional<Fraction> opt =
        computeOptimum(Direction::Up, con[constraintIndex]))
  return true;
return false;

[arjunp]

Addressed Siddharth's commments.

[ftynse]

We can't test for this because the Euclidean GCD algorithm is still correct and terminates with machine division (about efficiency, I am not sure).

You still can write a test computing the column echelon form for a given matrix that requires division by a negative value somewhere. The point of such a test is to (a) demonstrate that the result is correct and (b) provide enough coverage so that potential future changes in the column echelon computation affect a specific test rather than a much larger test (emptiness check) that can also be affected by dozens of other things and makes it ten times harder to find the root cause.

mlir/lib/Analysis/AffineStructures.cpp

1066

MLIR does not follow this pattern. If you prefer it had, please bring it up for general discussion rather than diverging locally.

[arjunp]

Added unit tests for makeTransformToColumnEchelon.

[ftynse]

Thank you!

[arjunp]

Great! I don't have commit access, can you land this? :)

mlir/include/mlir/Analysis/AffineStructures.h

 //===- AffineStructures.h - MLIR Affine Structures 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
 //
 //===----------------------------------------------------------------------===//
 //
 // Structures for affine/polyhedral analysis of ML functions.
 //
 //===----------------------------------------------------------------------===//
 
 #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"
 
 namespace mlir {
 
 class AffineCondition;
 class AffineForOp;
public:
 
   /// Runs the GCD test heuristic. If it proves inconclusive, falls back to
   /// generalized basis reduction if the set is bounded.
   ///
   /// Returns true if the set of constraints is found to have no solution,
   /// 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.
   ///
   /// Returns such a point if one exists, or an empty Optional otherwise.
   Optional<SmallVector<int64_t, 8>> findIntegerSample() const;
 
   /// Returns true if the given point satisfies the constraints, or false

mlir/include/mlir/Analysis/Presburger/Fraction.h

 inline Fraction operator-(Fraction x) { return Fraction(-x.num, x.den); }
 
 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; }
+
 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 Fraction operator*(Fraction x, Fraction y) {
   return Fraction(x.num * y.num, x.den * y.den);
 }
 
 } // namespace mlir
 
 #endif // MLIR_ANALYSIS_PRESBURGER_FRACTION_H

mlir/include/mlir/Analysis/Presburger/LinearTransform.h

+//===- 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_PRESBURGER_LINEARTRANSFORM_H
+#define MLIR_ANALYSIS_PRESBURGER_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.
+  static std::pair<unsigned, LinearTransform>
+  makeTransformToColumnEchelon(Matrix m);
+
+  // Post-multiply every constraint vector v in fac with this transform, say T,
+  // turning it into vT.
+  FlatAffineConstraints apply(const FlatAffineConstraints &fac);
+
+private:
+  // Post-multiply the given vector v with this transform, say T, returning vT.
+  SmallVector<int64_t, 8> postMultiplyVector(ArrayRef<int64_t> v);
+  explicit LinearTransform(Matrix oMatrix);
+
+  Matrix matrix;
+};
+
+} // namespace mlir
+#endif // MLIR_ANALYSIS_PRESBURGER_LINEARTRANSFORM_H

mlir/include/mlir/Analysis/Presburger/Matrix.h

public:
   unsigned getNumColumns() const;
 
   /// Get an ArrayRef corresponding to the specified row.
   ArrayRef<int64_t> getRow(unsigned row) const;
 
   /// 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.
   void resizeVertically(unsigned newNRows);
 
   /// Print the matrix.
   void print(raw_ostream &os) const;
   void dump() const;

mlir/include/mlir/Analysis/Presburger/Simplex.h

 //===----------------------------------------------------------------------===//
 
 #ifndef MLIR_ANALYSIS_PRESBURGER_SIMPLEX_H
 #define MLIR_ANALYSIS_PRESBURGER_SIMPLEX_H
 
 #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 {
 
 class GBRSimplex;
 
 /// This class implements a version of the Simplex and Generalized Basis
 /// Reduction algorithms, which can perform analysis of integer sets with affine

 /// column c, respectively.
 ///
 /// The unknowns in column position are together called the basis. Initially the
 /// basis is the set of variables -- in our example above, the initial basis is
 /// x, y.
 ///
 /// 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].
 ///
 /// The association of unknowns to rows and columns can be changed by a process
 /// called pivoting, where a row unknown and a column unknown exchange places
 /// and the remaining row variables' representation is changed accordingly
 /// by eliminating the old column unknown in favour of the new column unknown.
public:
 
   /// Rollback to a snapshot. This invalidates all later snapshots.
   void rollback(unsigned snapshot);
 
   /// Add all the constraints from the given FlatAffineConstraints.
   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<Fraction> 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<Fraction> computeOptimum(Direction direction,
                                     ArrayRef<int64_t> 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
   /// inequalities with successive indices.
   bool isMarkedRedundant(unsigned constraintIndex) const;
 
   /// Finds a subset of constraints that is redundant, i.e., such that
   /// the set of solutions does not change if these constraints are removed.
private:
   void swapRows(unsigned i, unsigned j);
 
   /// Restore the unknown to a non-negative sample value.
   ///
   /// Returns true if the unknown was successfully restored to a non-negative
   /// 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<Fraction> 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.
   void markRowRedundant(Unknown &u);
 
   /// Enum to denote operations that need to be undone during rollback.
   enum class UndoLogEntry {
     RemoveLastConstraint,

mlir/lib/Analysis/AffineStructures.cpp

 //===- AffineStructures.cpp - MLIR Affine Structures 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
 //
 //===----------------------------------------------------------------------===//
 //
 // Structures for affine/polyhedral analysis of affine dialect ops.
 //
 //===----------------------------------------------------------------------===//
 
 #include "mlir/Analysis/AffineStructures.h"
+#include "mlir/Analysis/Presburger/LinearTransform.h"
 #include "mlir/Analysis/Presburger/Simplex.h"
 #include "mlir/Dialect/Affine/IR/AffineOps.h"
 #include "mlir/Dialect/Affine/IR/AffineValueMap.h"
 #include "mlir/Dialect/StandardOps/IR/Ops.h"
 #include "mlir/IR/AffineExprVisitor.h"
 #include "mlir/IR/IntegerSet.h"
 #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"
 
 #define DEBUG_TYPE "affine-structures"
 
 using namespace mlir;
 using llvm::SmallDenseMap;
 using llvm::SmallDenseSet;
for (unsigned i = 0, e = getNumEqualities(); i < e; ++i) {
     int64_t v = std::abs(atEq(i, numCols - 1));
     if (gcd > 0 && (v % gcd != 0)) {
       return true;
     }
   }
   return false;
 }
 
-// First, try the GCD test heuristic.
-//
-// 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.
+// 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 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<unsigned, 8> 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.
+  Matrix dirs(boundedIneqs.size() + getNumEqualities(), getNumCols() - 1);
+
+  // Copy the bounded inequalities.
+  unsigned row = 0;
+  for (unsigned i : boundedIneqs) {
+    for (unsigned col = 0, e = getNumCols() - 1; col < e; ++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, e = getNumCols() - 1; col < e; ++col)
+      dirs(row, col) = atIneq(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) {
+  // unsigned unboundedDims = fac.getNumIds() - boundedDims;
+
+  // We iterate backwards so that whether we remove constraint i or not, the
+  // next constraint to be tested is always i - 1. i has to be signed to avoid
+  // overflow at zero.
+  for (int i = fac.getNumEqualities() - 1; i >= 0; i--)
+    if (eqInvolvesSuffixDims(fac, i, unboundedDims))
+      fac.removeEquality(i);
+  for (int i = fac.getNumInequalities() - 1; i >= 0; i--)
+    if (ineqInvolvesSuffixDims(fac, i, unboundedDims))
+      fac.removeInequality(i);
+}
+
+/// 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 an invertible 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 invertible, 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 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;
+  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 invertible as well, so this is a suitable
+  // transform to use in step 1 of our algorithm.
+  std::pair<unsigned, LinearTransform> result =
+      LinearTransform::makeTransformToColumnEchelon(std::move(m));
+  FlatAffineConstraints transformedSet = result.second.apply(*this);
+
+  unsigned numBoundedDims = result.first;
+  unsigned numUnboundedDims = getNumIds() - numBoundedDims;
+  removeConstraintsInvolvingSuffixDims(transformedSet, numUnboundedDims);
+  // Remove all the unbounded dimensions. i has to be signed to avoid overflow
+  // if boundedDims is zero.
+  for (int i = transformedSet.getNumIds() - 1; i >= int(numBoundedDims); i--)
+    transformedSet.removeId(numBoundedDims);
+
+  return !Simplex(transformedSet).findIntegerSample().hasValue();
 }
 
 Optional<SmallVector<int64_t, 8>>
 FlatAffineConstraints::findIntegerSample() const {
   return Simplex(*this).findIntegerSample();
 }
 
 /// Helper to evaluate an affine expression at a point.

mlir/lib/Analysis/Presburger/CMakeLists.txt

 add_mlir_library(MLIRPresburger
   Simplex.cpp
+  LinearTransform.cpp
   Matrix.cpp
-  )
+  )
 

mlir/lib/Analysis/Presburger/LinearTransform.cpp

+//===- 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/Presburger/LinearTransform.h"
+#include "mlir/Analysis/AffineStructures.h"
+
+namespace mlir {
+
+LinearTransform::LinearTransform(Matrix oMatrix) : matrix(std::move(oMatrix)) {}
+
+// Set M(row, targetCol) to its remainder on division by [0, M(row, sourceCol))
+// by subtracting from column targetCol an appropriate integer multiple of
+// sourceCol.
+// Apply the same column operation to otherMatrix. (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");
+  auto ratio = m(row, targetCol) / m(row, sourceCol);
+  m.addToColumn(sourceCol, targetCol, -ratio);
+  otherMatrix.addToColumn(sourceCol, targetCol, -ratio);
+}
+
+std::pair<unsigned, LinearTransform>
+LinearTransform::makeTransformToColumnEchelon(Matrix m) {
+  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` 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 on the same column, since this row is all zeros
+    // from column 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);
+    }
+
+    // Ensure the entry is non-negative.
+    if (m(row, echelonCol) < 0) {
+      m.negateColumn(echelonCol);
+      resultMatrix.negateColumn(echelonCol);
+    }
+
+    // Make all the other entries in `row` zero.
+    for (unsigned i = echelonCol + 1, e = m.getNumColumns(); i < e; ++i) {
+      // We 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 has value zero.
+
+      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)).
+      //
+      // 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) is not zero.
+      if (m(row, echelonCol) == 0) {
+        m.swapColumns(i, echelonCol);
+        resultMatrix.swapColumns(i, echelonCol);
+      }
+    }
+
+    ++echelonCol;
+  }
+
+  return {echelonCol, LinearTransform(std::move(resultMatrix))};
+}
+
+SmallVector<int64_t, 8>
+LinearTransform::postMultiplyVector(ArrayRef<int64_t> v) {
+  assert(v.size() == matrix.getNumRows() &&
+         "vector dimension should be matrix output dimension");
+
+  SmallVector<int64_t, 8> result;
+  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::apply(const FlatAffineConstraints &fac) {
+  FlatAffineConstraints result(fac.getNumDimIds());
+
+  for (unsigned i = 0, e = fac.getNumEqualities(); i < e; ++i) {
+    ArrayRef<int64_t> eq = fac.getEquality(i);
+
+    int64_t c = eq.back();
+
+    SmallVector<int64_t, 8> newEq = postMultiplyVector(eq.drop_back());
+    newEq.push_back(c);
+    result.addEquality(newEq);
+  }
+
+  for (unsigned i = 0, e = fac.getNumInequalities(); i < e; ++i) {
+    ArrayRef<int64_t> ineq = fac.getInequality(i);
+
+    int64_t c = ineq.back();
+
+    SmallVector<int64_t, 8> newIneq = postMultiplyVector(ineq.drop_back());
+    newIneq.push_back(c);
+    result.addInequality(newIneq);
+  }
+
+  return result;
+}
+
+} // namespace mlir
 No newline at end of file

mlir/lib/Analysis/Presburger/Matrix.cpp

 void Matrix::addToRow(unsigned sourceRow, unsigned targetRow, int64_t scale) {
   if (scale == 0)
     return;
   for (unsigned col = 0; col < nColumns; ++col)
     at(targetRow, col) += scale * at(sourceRow, col);
   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);
+  return;
+}
+
+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)
       os << at(row, column) << ' ';
     os << '\n';
   }
 }
 
 void Matrix::dump() const { print(llvm::errs()); }
 
 } // namespace mlir

mlir/lib/Analysis/Presburger/Simplex.cpp

 //===- Simplex.cpp - MLIR Simplex Class -----------------------------------===//
 //
 // Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
 // See https://llvm.org/LICENSE.txt for license information.
 // SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
 //
 //===----------------------------------------------------------------------===//
 
 #include "mlir/Analysis/Presburger/Simplex.h"
 #include "mlir/Analysis/Presburger/Matrix.h"
 #include "mlir/Support/MathExtras.h"
+#include "llvm/ADT/Optional.h"
 
 namespace mlir {
 using Direction = Simplex::Direction;
 
 const int nullIndex = std::numeric_limits<int>::max();
 
 /// Construct a Simplex object with `nVar` variables.
 Simplex::Simplex(unsigned nVar)
Optional<Fraction> Simplex::computeRowOptimum(Direction direction,
   // denominator for this row.
   return Fraction(tableau(row, 1), tableau(row, 0));
 }
 
 /// Compute the optimum of the specified expression in the specified direction,
 /// or None if it is unbounded.
 Optional<Fraction> Simplex::computeOptimum(Direction direction,
                                            ArrayRef<int64_t> coeffs) {
-  assert(!empty && "Tableau should not be empty");
+  assert(!empty && "Simplex should not be empty");
 
   unsigned snapshot = getSnapshot();
   unsigned conIndex = addRow(coeffs);
   unsigned row = con[conIndex].pos;
   Optional<Fraction> optimum = computeRowOptimum(direction, row);
   rollback(snapshot);
   return optimum;
 }
 
+Optional<Fraction> Simplex::computeOptimum(Direction direction, Unknown &u) {
+  assert(!empty && "Simplex should not be empty!");
+  if (u.orientation == Orientation::Column) {
+    unsigned column = u.pos;
+    Optional<unsigned> 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<Fraction> 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.
+  if (Optional<Fraction> opt =
+          computeOptimum(Direction::Up, con[constraintIndex]))
+    return true;
+  return false;
+}
+
 /// Redundant constraints are those that are in row orientation and lie in
 /// rows 0 to nRedundant - 1.
 bool Simplex::isMarkedRedundant(unsigned constraintIndex) const {
   const Unknown &u = con[constraintIndex];
   return u.orientation == Orientation::Row && u.pos < nRedundant;
 }
 
 /// Mark the specified row redundant.

mlir/unittests/Analysis/AffineStructuresTest.cpp

 
 #include <gmock/gmock.h>
 #include <gtest/gtest.h>
 
 #include <numeric>
 
 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<SmallVector<int64_t, 8>> maybeSample = fac.findIntegerSample();
-  if (!hasValue) {
+/// 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<SmallVector<int64_t, 8>> 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));
     }
+    break;
+  case TestFunction::Empty:
+    EXPECT_EQ(!hasSample, fac.isIntegerEmpty());
+    break;
+  }
 }
 
 /// Construct a FlatAffineConstraints from a set of inequality and
 /// equality constraints.
 static FlatAffineConstraints
 makeFACFromConstraints(unsigned dims, ArrayRef<SmallVector<int64_t, 4>> ineqs,
                        ArrayRef<SmallVector<int64_t, 4>> eqs) {
   FlatAffineConstraints fac(ineqs.size(), eqs.size(), dims + 1, dims);
   for (const auto &eq : eqs)
     fac.addEquality(eq);
   for (const auto &ineq : ineqs)
     fac.addInequality(ineq);
   return fac;
 }
 
 /// Check sampling for all the permutations of the dimensions for the given
 /// constraint set. Since the GBR algorithm progresses dimension-wise, different
 /// 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<SmallVector<int64_t, 4>> ineqs,
-                                    ArrayRef<SmallVector<int64_t, 4>> eqs) {
+                                    ArrayRef<SmallVector<int64_t, 4>> eqs,
+                                    TestFunction fn = TestFunction::Sample) {
   SmallVector<unsigned, 4> perm(nDim);
   std::iota(perm.begin(), perm.end(), 0);
   auto permute = [&perm](ArrayRef<int64_t> coeffs) {
     SmallVector<int64_t, 4> permuted;
     for (unsigned id : perm)
       permuted.push_back(coeffs[id]);
     permuted.push_back(coeffs.back());
     return permuted;
   };
   do {
     SmallVector<SmallVector<int64_t, 4>, 4> permutedIneqs, permutedEqs;
     for (const auto &ineq : ineqs)
       permutedIneqs.push_back(permute(ineq));
     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()));
 }
 
 TEST(FlatAffineConstraintsTest, FindSampleTest) {
   // Bounded sets with only inequalities.
 
   // 0 <= 7x <= 5
   checkSample(true, makeFACFromConstraints(1, {{7, 0}, {-7, 5}}, {}));

 TEST(FlatAffineConstraintsTest, IsIntegerEmptyTest) {
   // 1 <= 5x and 5x <= 4 (no solution).
   EXPECT_TRUE(
       makeFACFromConstraints(1, {{5, -1}, {-5, 4}}, {}).isIntegerEmpty());
   // 1 <= 5x and 5x <= 9 (solution: x = 1).
   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 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);
+
   // FlatAffineConstraints::isEmpty() does not detect the following sets to be
   // empty.
 
   // 3x + 7y = 1 and 0 <= x, y <= 10.
   // Since x and y are non-negative, 3x + 7y can never be 1.
   EXPECT_TRUE(
       makeFACFromConstraints(
           2, {{1, 0, 0}, {-1, 0, 10}, {0, 1, 0}, {0, -1, 10}}, {{3, 7, -1}})