[MLIR] Move presburger math from FlatAffineConstraints to Presburger directory

[Groverkss]

Summary

This patch moves remaining math functionality in FlatAffineConstraints to
IntegerPolyhedron. This patch also moves LinearTransform to Presburger/
directory and renames FlatAffineConstraints to IntegerPolyhedron in
Presburger/Simplex and Presburger/Utils.

This patch is a continuation of https://reviews.llvm.org/D114674.

[Groverkss]

The only non-trivial change in this patch that would require you to look at the implementation is IntegerPolyhedron::print. FlatAffineValueConstraints prints information about SSA Values, this information is now printed by a virtual function IntegerPolyhedron::printSpace which FlatAffineConstraints overrides to include the SSA Value information.

The rest of the patch is only moving and renamed to the Presburger directory.

[Groverkss]

• Fix clang-format error

[arjunp]

I would suggest just doing the move of functions from FAC to IP in this patch, and doing the other things like using IP in PresbugerSet and LinearTransform and moving LinearTransform in different patches.

Also, in the summary I would change "renames FlatAffineConstraints to IntegerPolyhedron" to "uses "IntegerPolyhedron instead of FlatAffineConstraints".

mlir/include/mlir/Analysis/AffineStructures.h

115

Why is this needed exactly? Might be useful to write a comment.

131

nit: seems to be moved unnecessarily

mlir/include/mlir/Analysis/PresburgerSet.h

120

mlir/lib/Analysis/AffineStructures.cpp

1344–1347

Can't you dyn_cast the IntegerPolyhedron directly to a FlatAffineValueConstraints?

1351

Isn't this already a FlatAffineConstraints *?

Is the static cast to other needed here?

1356–1362

Again is otherFAC needed to be used here?

mlir/lib/Analysis/Presburger/Simplex.cpp

507

I think variable names like fac should be changed throughout.

mlir/lib/Analysis/PresburgerSet.cpp

29

424

[Groverkss]

Thank you for the review. I'll split this patch up into smaller patches to make it easier to review and take the comments in consideration.

mlir/include/mlir/Analysis/AffineStructures.h

public:
 
   virtual ~FlatAffineConstraints() = default;
 
   /// Return the kind of this FlatAffineConstraints.
   virtual Kind getKind() const { return Kind::FlatAffineConstraints; }
 
   static bool classof(const FlatAffineConstraints *cst) { return true; }
 
-  /// Checks for emptiness by performing variable elimination on all
-  /// identifiers, running the GCD test on each equality constraint, and
-  /// checking for invalid constraints. Returns true if the GCD test fails for
-  /// any equality, or if any invalid constraints are discovered on any row.
-  /// Returns false otherwise.
-  bool isEmpty() const;
-
-  /// Runs the GCD test on all equality constraints. Returns true if this test
-  /// fails on any equality. Returns false otherwise.
-  /// This test can be used to disprove the existence of a solution. If it
-  /// returns true, no integer solution to the equality constraints can exist.
-  bool isEmptyByGCDTest() const;
-
-  /// Returns true if the set of constraints is found to have no solution,
-  /// false if a solution exists. Uses the same algorithm as
-  /// `findIntegerSample`.
-  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 an integer sample point satisfying the constraints using a
-  /// branch and bound algorithm with generalized basis reduction, with some
-  /// additional processing using Simplex for unbounded sets.
-  ///
-  /// Returns an integer sample 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
-  /// otherwise.
-  bool containsPoint(ArrayRef<int64_t> point) const;
-
-  /// Find pairs of inequalities identified by their position indices, using
-  /// which an explicit representation for each local variable can be computed.
-  /// The pairs are stored as indices of upperbound, lowerbound inequalities. If
-  /// no such pair can be found, it is stored as llvm::None.
-  ///
-  /// The dividends of the explicit representations are stored in `dividends`
-  /// and the denominators in `denominators`. If no explicit representation
-  /// could be found for the `i^th` local identifier, `denominators[i]` is set
-  /// to 0.
-  void getLocalReprs(
-      std::vector<SmallVector<int64_t, 8>> &dividends,
-      SmallVector<unsigned, 4> &denominators,
-      std::vector<llvm::Optional<std::pair<unsigned, unsigned>>> &repr) const;
-  void getLocalReprs(
-      std::vector<llvm::Optional<std::pair<unsigned, unsigned>>> &repr) const;
-  void getLocalReprs(std::vector<SmallVector<int64_t, 8>> &dividends,
-                     SmallVector<unsigned, 4> &denominators) const;
-
   // Clones this object.
   std::unique_ptr<FlatAffineConstraints> clone() const;
 
-  /// The type of bound: equal, lower bound or upper bound.
-  enum BoundType { EQ, LB, UB };
-
   /// Adds a bound for the identifier at the specified position with constraints
   /// being drawn from the specified bound map. In case of an EQ bound, the
   /// bound map is expected to have exactly one result. In case of a LB/UB, the
   /// bound map may have more than one result, for each of which an inequality
   /// is added.
   /// Note: The dimensions/symbols of this FlatAffineConstraints must match the
   /// dimensions/symbols of the affine map.
   LogicalResult addBound(BoundType type, unsigned pos, AffineMap boundMap);
 
-  /// Adds a constant bound for the specified identifier.
-  void addBound(BoundType type, unsigned pos, int64_t value);
-
-  /// Adds a constant bound for the specified expression.
-  void addBound(BoundType type, ArrayRef<int64_t> expr, int64_t value);
-
-  /// Returns the constraint system as an integer set. Returns a null integer
-  /// set if the system has no constraints, or if an integer set couldn't be
-  /// constructed as a result of a local variable's explicit representation not
-  /// being known and such a local variable appearing in any of the constraints.
-  IntegerSet getAsIntegerSet(MLIRContext *context) const;
+  using IntegerPolyhedron::addBound;
 
   /// Computes the lower and upper bounds of the first `num` dimensional
   /// identifiers (starting at `offset`) as an affine map of the remaining
   /// identifiers (dimensional and symbolic). This method is able to detect
   /// identifiers as floordiv's and mod's of affine expressions of other
   /// identifiers with respect to (positive) constants. Sets bound map to a
   /// null AffineMap if such a bound can't be found (or yet unimplemented).
   void getSliceBounds(unsigned offset, unsigned num, MLIRContext *context,
                       SmallVectorImpl<AffineMap> *lbMaps,
                       SmallVectorImpl<AffineMap> *ubMaps);
 
-  /// Adds a new local identifier as the floordiv of an affine function of other
-  /// identifiers, the coefficients of which are provided in `dividend` and with
-  /// respect to a positive constant `divisor`. Two constraints are added to the
-  /// system to capture equivalence with the floordiv:
-  /// q = dividend floordiv c    <=>   c*q <= dividend <= c*q + c - 1.
-  void addLocalFloorDiv(ArrayRef<int64_t> dividend, int64_t divisor);
+  /// Returns the constraint system as an integer set. Returns a null integer
+  /// set if the system has no constraints, or if an integer set couldn't be
+  /// constructed as a result of a local variable's explicit representation not
+  /// being known and such a local variable appearing in any of the constraints.
+  IntegerSet getAsIntegerSet(MLIRContext *context) const;
 
   /// Composes an affine map whose dimensions and symbols match one to one with
   /// the dimensions and symbols of this FlatAffineConstraints. The results of
   /// the map `other` are added as the leading dimensions of this constraint
   /// system. Returns failure if `other` is a semi-affine map.
   LogicalResult composeMatchingMap(AffineMap other);
 
-  /// Projects out (aka eliminates) `num` identifiers starting at position
-  /// `pos`. The resulting constraint system is the shadow along the dimensions
-  /// that still exist. This method may not always be integer exact.
-  // TODO: deal with integer exactness when necessary - can return a value to
-  // mark exactness for example.
-  void projectOut(unsigned pos, unsigned num);
-  inline void projectOut(unsigned pos) { return projectOut(pos, 1); }
-
-  /// Sets the `values.size()` identifiers starting at `po`s to the specified
-  /// values and removes them.
-  void setAndEliminate(unsigned pos, ArrayRef<int64_t> values);
-
-  /// Changes the partition between dimensions and symbols. Depending on the new
-  /// symbol count, either a chunk of trailing dimensional identifiers becomes
-  /// symbols, or some of the leading symbols become dimensions.
-  void setDimSymbolSeparation(unsigned newSymbolCount);
-
-  /// Tries to fold the specified identifier to a constant using a trivial
-  /// equality detection; if successful, the constant is substituted for the
-  /// identifier everywhere in the constraint system and then removed from the
-  /// system.
-  LogicalResult constantFoldId(unsigned pos);
-
-  /// This method calls `constantFoldId` for the specified range of identifiers,
-  /// `num` identifiers starting at position `pos`.
-  void constantFoldIdRange(unsigned pos, unsigned num);
-
-  /// Updates the constraints to be the smallest bounding (enclosing) box that
-  /// contains the points of `this` set and that of `other`, with the symbols
-  /// being treated specially. For each of the dimensions, the min of the lower
-  /// bounds (symbolic) and the max of the upper bounds (symbolic) is computed
-  /// to determine such a bounding box. `other` is expected to have the same
-  /// dimensional identifiers as this constraint system (in the same order).
-  ///
-  /// E.g.:
-  /// 1) this   = {0 <= d0 <= 127},
-  ///    other  = {16 <= d0 <= 192},
-  ///    output = {0 <= d0 <= 192}
-  /// 2) this   = {s0 + 5 <= d0 <= s0 + 20},
-  ///    other  = {s0 + 1 <= d0 <= s0 + 9},
-  ///    output = {s0 + 1 <= d0 <= s0 + 20}
-  /// 3) this   = {0 <= d0 <= 5, 1 <= d1 <= 9}
-  ///    other  = {2 <= d0 <= 6, 5 <= d1 <= 15},
-  ///    output = {0 <= d0 <= 6, 1 <= d1 <= 15}
-  LogicalResult unionBoundingBox(const FlatAffineConstraints &other);
-
   /// Replaces the contents of this FlatAffineConstraints with `other`.
-  virtual void clearAndCopyFrom(const FlatAffineConstraints &other);
-
-  /// Returns the smallest known constant bound for the extent of the specified
-  /// identifier (pos^th), i.e., the smallest known constant that is greater
-  /// than or equal to 'exclusive upper bound' - 'lower bound' of the
-  /// identifier. This constant bound is guaranteed to be non-negative. Returns
-  /// None if it's not a constant. This method employs trivial (low complexity /
-  /// cost) checks and detection. Symbolic identifiers are treated specially,
-  /// i.e., it looks for constant differences between affine expressions
-  /// involving only the symbolic identifiers. `lb` and `ub` (along with the
-  /// `boundFloorDivisor`) are set to represent the lower and upper bound
-  /// associated with the constant difference: `lb`, `ub` have the coefficients,
-  /// and `boundFloorDivisor`, their divisor. `minLbPos` and `minUbPos` if
-  /// non-null are set to the position of the constant lower bound and upper
-  /// bound respectively (to the same if they are from an equality). Ex: if the
-  /// lower bound is [(s0 + s2 - 1) floordiv 32] for a system with three
-  /// symbolic identifiers, *lb = [1, 0, 1], lbDivisor = 32. See comments at
-  /// function definition for examples.
-  Optional<int64_t> getConstantBoundOnDimSize(
-      unsigned pos, SmallVectorImpl<int64_t> *lb = nullptr,
-      int64_t *boundFloorDivisor = nullptr,
-      SmallVectorImpl<int64_t> *ub = nullptr, unsigned *minLbPos = nullptr,
-      unsigned *minUbPos = nullptr) const;
-
-  /// Returns the constant bound for the pos^th identifier if there is one;
-  /// None otherwise.
-  // TODO: Support EQ bounds.
-  Optional<int64_t> getConstantBound(BoundType type, unsigned pos) const;
+  virtual void clearAndCopyFrom(const IntegerPolyhedron &other) override;
 
   /// Gets the lower and upper bound of the `offset` + `pos`th identifier
   /// treating [0, offset) U [offset + num, symStartPos) as dimensions and
   /// [symStartPos, getNumDimAndSymbolIds) as symbols, and `pos` lies in
   /// [0, num). The multi-dimensional maps in the returned pair represent the
   /// max and min of potentially multiple affine expressions. The upper bound is
   /// exclusive. `localExprs` holds pre-computed AffineExpr's for all local
   /// identifiers in the system.
   std::pair<AffineMap, AffineMap>
   getLowerAndUpperBound(unsigned pos, unsigned offset, unsigned num,
                         unsigned symStartPos, ArrayRef<AffineExpr> localExprs,
                         MLIRContext *context) const;
 
-  /// Removes constraints that are independent of (i.e., do not have a
-  /// coefficient) identifiers in the range [pos, pos + num).
-  void removeIndependentConstraints(unsigned pos, unsigned num);
-
-  /// Returns true if the set can be trivially detected as being
-  /// hyper-rectangular on the specified contiguous set of identifiers.
-  bool isHyperRectangular(unsigned pos, unsigned num) const;
-
-  /// Removes duplicate constraints, trivially true constraints, and constraints
-  /// that can be detected as redundant as a result of differing only in their
-  /// constant term part. A constraint of the form <non-negative constant> >= 0
-  /// is considered trivially true. This method is a linear time method on the
-  /// constraints, does a single scan, and updates in place. It also normalizes
-  /// constraints by their GCD and performs GCD tightening on inequalities.
-  void removeTrivialRedundancy();
-
-  /// A more expensive check than `removeTrivialRedundancy` to detect redundant
-  /// inequalities.
-  void removeRedundantInequalities();
-
-  /// Removes redundant constraints using Simplex. Although the algorithm can
-  /// theoretically take exponential time in the worst case (rare), it is known
-  /// to perform much better in the average case. If V is the number of vertices
-  /// in the polytope and C is the number of constraints, the algorithm takes
-  /// O(VC) time.
-  void removeRedundantConstraints();
-
-  /// Converts identifiers in the column range [idStart, idLimit) to local
-  /// variables.
-  void convertDimToLocal(unsigned dimStart, unsigned dimLimit);
-
-  /// Adds additional local ids to the sets such that they both have the union
-  /// of the local ids in each set, without changing the set of points that
-  /// lie in `this` and `other`. The ordering of the local ids in the
-  /// sets may also be changed. After merging, if the `i^th` local variable in
-  /// one set has a known division representation, then the `i^th` local
-  /// variable in the other set either has the same division representation or
-  /// no known division representation.
-  ///
-  /// The number of dimensions and symbol ids in `this` and `other` should
-  /// match.
-  void mergeLocalIds(FlatAffineConstraints &other);
-
-  void print(raw_ostream &os) const;
-  void dump() const;
-
 protected:
-  /// Returns false if the fields corresponding to various identifier counts, or
-  /// equality/inequality buffer sizes aren't consistent; true otherwise. This
-  /// is meant to be used within an assert internally.
-  virtual bool hasConsistentState() const;
-
-  /// Checks all rows of equality/inequality constraints for trivial
-  /// contradictions (for example: 1 == 0, 0 >= 1), which may have surfaced
-  /// after elimination. Returns true if an invalid constraint is found;
-  /// false otherwise.
-  bool hasInvalidConstraint() const;
-
-  /// Returns the constant lower bound bound if isLower is true, and the upper
-  /// bound if isLower is false.
-  template <bool isLower>
-  Optional<int64_t> computeConstantLowerOrUpperBound(unsigned pos);
-
   /// Given an affine map that is aligned with this constraint system:
   /// * Flatten the map.
   /// * Add newly introduced local columns at the beginning of this constraint
   ///   system (local column pos 0).
   /// * Add equalities that define the new local columns to this constraint
   ///   system.
   /// * Return the flattened expressions via `flattenedExprs`.
   ///
   /// Note: This is a shared helper function of `addLowerOrUpperBound` and
   ///       `composeMatchingMap`.
   LogicalResult flattenAlignedMapAndMergeLocals(
       AffineMap map, std::vector<SmallVector<int64_t, 8>> *flattenedExprs);
 
-  /// Eliminates a single identifier at `position` from equality and inequality
-  /// constraints. Returns `success` if the identifier was eliminated, and
-  /// `failure` otherwise.
-  inline LogicalResult gaussianEliminateId(unsigned position) {
-    return success(gaussianEliminateIds(position, position + 1) == 1);
-  }
-
-  /// Removes local variables using equalities. Each equality is checked if it
-  /// can be reduced to the form: `e = affine-expr`, where `e` is a local
-  /// variable and `affine-expr` is an affine expression not containing `e`.
-  /// If an equality satisfies this form, the local variable is replaced in
-  /// each constraint and then removed. The equality used to replace this local
-  /// variable is also removed.
-  void removeRedundantLocalVars();
-
-  /// Eliminates identifiers from equality and inequality constraints
-  /// in column range [posStart, posLimit).
-  /// Returns the number of variables eliminated.
-  unsigned gaussianEliminateIds(unsigned posStart, unsigned posLimit);
-
-  /// Eliminates the identifier at the specified position using Fourier-Motzkin
-  /// variable elimination, but uses Gaussian elimination if there is an
-  /// equality involving that identifier. If the result of the elimination is
-  /// integer exact, `*isResultIntegerExact` is set to true. If `darkShadow` is
-  /// set to true, a potential under approximation (subset) of the rational
-  /// shadow / exact integer shadow is computed.
-  // See implementation comments for more details.
-  virtual void fourierMotzkinEliminate(unsigned pos, bool darkShadow = false,
-                                       bool *isResultIntegerExact = nullptr);
-
-  /// Tightens inequalities given that we are dealing with integer spaces. This
-  /// is similar to the GCD test but applied to inequalities. The constant term
-  /// can be reduced to the preceding multiple of the GCD of the coefficients,
-  /// i.e.,
-  ///  64*i - 100 >= 0  =>  64*i - 128 >= 0 (since 'i' is an integer). This is a
-  /// fast method (linear in the number of coefficients).
-  void gcdTightenInequalities();
-
-  /// Normalized each constraints by the GCD of its coefficients.
-  void normalizeConstraintsByGCD();
-
-  /// Searches for a constraint with a non-zero coefficient at `colIdx` in
-  /// equality (isEq=true) or inequality (isEq=false) constraints.
-  /// Returns true and sets row found in search in `rowIdx`, false otherwise.
-  bool findConstraintWithNonZeroAt(unsigned colIdx, bool isEq,
-                                   unsigned *rowIdx) const;
-
-  /// Returns true if the pos^th column is all zero for both inequalities and
-  /// equalities.
-  bool isColZero(unsigned pos) const;
-
-  /// A parameter that controls detection of an unrealistic number of
-  /// constraints. If the number of constraints is this many times the number of
-  /// variables, we consider such a system out of line with the intended use
-  /// case of FlatAffineConstraints.
-  // The rationale for 32 is that in the typical simplest of cases, an
-  // identifier is expected to have one lower bound and one upper bound
-  // constraint. With a level of tiling or a connection to another identifier
-  // through a div or mod, an extra pair of bounds gets added. As a limit, we
-  // don't expect an identifier to have more than 32 lower/upper/equality
-  // constraints. This is conservatively set low and can be raised if needed.
-  constexpr static unsigned kExplosionFactor = 32;
+  /// Prints the number of constraints, dimensions, symbols, and locals. Also
+  /// prints for each identifier, if there is a SSA value attached to it.
+  void printSpace(raw_ostream &os) const override;
 };
 
 /// An extension of FlatAffineConstraints in which dimensions and symbols can
 /// optionally be associated with an SSA value.
 class FlatAffineValueConstraints : public FlatAffineConstraints {
 public:
   /// Constructs a constraint system reserving memory for the specified number
   /// of constraints and identifiers.
void mergeAndAlignIdsWithOther(unsigned offset,
                                  FlatAffineValueConstraints *other);
 
   /// Returns true if this constraint system and `other` are in the same
   /// space, i.e., if they are associated with the same set of identifiers,
   /// appearing in the same order. Returns false otherwise.
   bool areIdsAlignedWithOther(const FlatAffineValueConstraints &other);
 
   /// Replaces the contents of this FlatAffineValueConstraints with `other`.
-  void clearAndCopyFrom(const FlatAffineConstraints &other) override;
+  void clearAndCopyFrom(const IntegerPolyhedron &other) override;
 
   /// Returns the Value associated with the pos^th identifier. Asserts if
   /// no Value identifier was associated.
   inline Value getValue(unsigned pos) const {
     assert(hasValue(pos) && "identifier's Value not set");
     return values[pos].getValue();
   }
 

mlir/include/mlir/Analysis/Presburger/IntegerPolyhedron.h

 //===- IntegerPolyhedron.h - MLIR IntegerPolyhedron 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
 //
 //===----------------------------------------------------------------------===//
 //
 // A class to represent an integer polyhedron.
 //
 //===----------------------------------------------------------------------===//
 
 #ifndef MLIR_ANALYSIS_PRESBURGER_INTEGERPOLYHEDRON_H
 #define MLIR_ANALYSIS_PRESBURGER_INTEGERPOLYHEDRON_H
 
 #include "mlir/Analysis/Presburger/Matrix.h"
+#include "mlir/Support/LogicalResult.h"
 
 namespace mlir {
 
 /// An integer polyhedron is the set of solutions to a list of affine
 /// constraints over n integer-valued variables/identifiers. Affine constraints
 /// can be inequalities or equalities in the form:
 ///
 /// Inequality: c_0*x_0 + c_1*x_1 + .... + c_{n-1}*x_{n-1} + c_n >= 0
public:
   void removeInequalityRange(unsigned start, unsigned end);
 
   /// Swap the posA^th identifier with the posB^th identifier.
   virtual void swapId(unsigned posA, unsigned posB);
 
   /// Removes all equalities and inequalities.
   void clearConstraints();
 
+  /// Return a system with no constraints, i.e., one which is satisfied by all
+  /// points.
+  static IntegerPolyhedron getUniverse(unsigned numDims = 0,
+                                       unsigned numSymbols = 0) {
+    return IntegerPolyhedron(numDims, numSymbols);
+  }
+
+  /// Checks for emptiness by performing variable elimination on all
+  /// identifiers, running the GCD test on each equality constraint, and
+  /// checking for invalid constraints. Returns true if the GCD test fails for
+  /// any equality, or if any invalid constraints are discovered on any row.
+  /// Returns false otherwise.
+  bool isEmpty() const;
+
+  /// Runs the GCD test on all equality constraints. Returns true if this test
+  /// fails on any equality. Returns false otherwise.
+  /// This test can be used to disprove the existence of a solution. If it
+  /// returns true, no integer solution to the equality constraints can exist.
+  bool isEmptyByGCDTest() const;
+
+  /// Returns true if the set of constraints is found to have no solution,
+  /// false if a solution exists. Uses the same algorithm as
+  /// `findIntegerSample`.
+  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 an integer sample point satisfying the constraints using a
+  /// branch and bound algorithm with generalized basis reduction, with some
+  /// additional processing using Simplex for unbounded sets.
+  ///
+  /// Returns an integer sample 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
+  /// otherwise.
+  bool containsPoint(ArrayRef<int64_t> point) const;
+
+  /// Find pairs of inequalities identified by their position indices, using
+  /// which an explicit representation for each local variable can be computed.
+  /// The pairs are stored as indices of upperbound, lowerbound inequalities. If
+  /// no such pair can be found, it is stored as llvm::None.
+  ///
+  /// The dividends of the explicit representations are stored in `dividends`
+  /// and the denominators in `denominators`. If no explicit representation
+  /// could be found for the `i^th` local identifier, `denominators[i]` is set
+  /// to 0.
+  void getLocalReprs(
+      std::vector<SmallVector<int64_t, 8>> &dividends,
+      SmallVector<unsigned, 4> &denominators,
+      std::vector<llvm::Optional<std::pair<unsigned, unsigned>>> &repr) const;
+  void getLocalReprs(
+      std::vector<llvm::Optional<std::pair<unsigned, unsigned>>> &repr) const;
+  void getLocalReprs(std::vector<SmallVector<int64_t, 8>> &dividends,
+                     SmallVector<unsigned, 4> &denominators) const;
+
+  /// The type of bound: equal, lower bound or upper bound.
+  enum BoundType { EQ, LB, UB };
+
+  /// Adds a constant bound for the specified identifier.
+  void addBound(BoundType type, unsigned pos, int64_t value);
+
+  /// Adds a constant bound for the specified expression.
+  void addBound(BoundType type, ArrayRef<int64_t> expr, int64_t value);
+
+  /// Adds a new local identifier as the floordiv of an affine function of other
+  /// identifiers, the coefficients of which are provided in `dividend` and with
+  /// respect to a positive constant `divisor`. Two constraints are added to the
+  /// system to capture equivalence with the floordiv:
+  /// q = dividend floordiv c    <=>   c*q <= dividend <= c*q + c - 1.
+  void addLocalFloorDiv(ArrayRef<int64_t> dividend, int64_t divisor);
+
+  /// Projects out (aka eliminates) `num` identifiers starting at position
+  /// `pos`. The resulting constraint system is the shadow along the dimensions
+  /// that still exist. This method may not always be integer exact.
+  // TODO: deal with integer exactness when necessary - can return a value to
+  // mark exactness for example.
+  void projectOut(unsigned pos, unsigned num);
+  inline void projectOut(unsigned pos) { return projectOut(pos, 1); }
+
+  /// Sets the `values.size()` identifiers starting at `po`s to the specified
+  /// values and removes them.
+  void setAndEliminate(unsigned pos, ArrayRef<int64_t> values);
+
+  /// Changes the partition between dimensions and symbols. Depending on the new
+  /// symbol count, either a chunk of trailing dimensional identifiers becomes
+  /// symbols, or some of the leading symbols become dimensions.
+  void setDimSymbolSeparation(unsigned newSymbolCount);
+
+  /// Tries to fold the specified identifier to a constant using a trivial
+  /// equality detection; if successful, the constant is substituted for the
+  /// identifier everywhere in the constraint system and then removed from the
+  /// system.
+  LogicalResult constantFoldId(unsigned pos);
+
+  /// This method calls `constantFoldId` for the specified range of identifiers,
+  /// `num` identifiers starting at position `pos`.
+  void constantFoldIdRange(unsigned pos, unsigned num);
+
+  /// Updates the constraints to be the smallest bounding (enclosing) box that
+  /// contains the points of `this` set and that of `other`, with the symbols
+  /// being treated specially. For each of the dimensions, the min of the lower
+  /// bounds (symbolic) and the max of the upper bounds (symbolic) is computed
+  /// to determine such a bounding box. `other` is expected to have the same
+  /// dimensional identifiers as this constraint system (in the same order).
+  ///
+  /// E.g.:
+  /// 1) this   = {0 <= d0 <= 127},
+  ///    other  = {16 <= d0 <= 192},
+  ///    output = {0 <= d0 <= 192}
+  /// 2) this   = {s0 + 5 <= d0 <= s0 + 20},
+  ///    other  = {s0 + 1 <= d0 <= s0 + 9},
+  ///    output = {s0 + 1 <= d0 <= s0 + 20}
+  /// 3) this   = {0 <= d0 <= 5, 1 <= d1 <= 9}
+  ///    other  = {2 <= d0 <= 6, 5 <= d1 <= 15},
+  ///    output = {0 <= d0 <= 6, 1 <= d1 <= 15}
+  LogicalResult unionBoundingBox(const IntegerPolyhedron &other);
+
+  /// Returns the smallest known constant bound for the extent of the specified
+  /// identifier (pos^th), i.e., the smallest known constant that is greater
+  /// than or equal to 'exclusive upper bound' - 'lower bound' of the
+  /// identifier. This constant bound is guaranteed to be non-negative. Returns
+  /// None if it's not a constant. This method employs trivial (low complexity /
+  /// cost) checks and detection. Symbolic identifiers are treated specially,
+  /// i.e., it looks for constant differences between affine expressions
+  /// involving only the symbolic identifiers. `lb` and `ub` (along with the
+  /// `boundFloorDivisor`) are set to represent the lower and upper bound
+  /// associated with the constant difference: `lb`, `ub` have the coefficients,
+  /// and `boundFloorDivisor`, their divisor. `minLbPos` and `minUbPos` if
+  /// non-null are set to the position of the constant lower bound and upper
+  /// bound respectively (to the same if they are from an equality). Ex: if the
+  /// lower bound is [(s0 + s2 - 1) floordiv 32] for a system with three
+  /// symbolic identifiers, *lb = [1, 0, 1], lbDivisor = 32. See comments at
+  /// function definition for examples.
+  Optional<int64_t> getConstantBoundOnDimSize(
+      unsigned pos, SmallVectorImpl<int64_t> *lb = nullptr,
+      int64_t *boundFloorDivisor = nullptr,
+      SmallVectorImpl<int64_t> *ub = nullptr, unsigned *minLbPos = nullptr,
+      unsigned *minUbPos = nullptr) const;
+
+  /// Returns the constant bound for the pos^th identifier if there is one;
+  /// None otherwise.
+  // TODO: Support EQ bounds.
+  Optional<int64_t> getConstantBound(BoundType type, unsigned pos) const;
+
+  /// Replaces the contents of this IntegerPolyhedron with `other`.
+  virtual void clearAndCopyFrom(const IntegerPolyhedron &other);
+
   /// Gather positions of all lower and upper bounds of the identifier at `pos`,
   /// and optionally any equalities on it. In addition, the bounds are to be
   /// independent of identifiers in position range [`offset`, `offset` + `num`).
   void
   getLowerAndUpperBoundIndices(unsigned pos,
                                SmallVectorImpl<unsigned> *lbIndices,
                                SmallVectorImpl<unsigned> *ubIndices,
                                SmallVectorImpl<unsigned> *eqIndices = nullptr,
                                unsigned offset = 0, unsigned num = 0) const;
 
+  /// Removes constraints that are independent of (i.e., do not have a
+  /// coefficient) identifiers in the range [pos, pos + num).
+  void removeIndependentConstraints(unsigned pos, unsigned num);
+
+  /// Returns true if the set can be trivially detected as being
+  /// hyper-rectangular on the specified contiguous set of identifiers.
+  bool isHyperRectangular(unsigned pos, unsigned num) const;
+
+  /// Removes duplicate constraints, trivially true constraints, and constraints
+  /// that can be detected as redundant as a result of differing only in their
+  /// constant term part. A constraint of the form <non-negative constant> >= 0
+  /// is considered trivially true. This method is a linear time method on the
+  /// constraints, does a single scan, and updates in place. It also normalizes
+  /// constraints by their GCD and performs GCD tightening on inequalities.
+  void removeTrivialRedundancy();
+
+  /// A more expensive check than `removeTrivialRedundancy` to detect redundant
+  /// inequalities.
+  void removeRedundantInequalities();
+
+  /// Removes redundant constraints using Simplex. Although the algorithm can
+  /// theoretically take exponential time in the worst case (rare), it is known
+  /// to perform much better in the average case. If V is the number of vertices
+  /// in the polytope and C is the number of constraints, the algorithm takes
+  /// O(VC) time.
+  void removeRedundantConstraints();
+
+  /// Converts identifiers in the column range [idStart, idLimit) to local
+  /// variables.
+  void convertDimToLocal(unsigned dimStart, unsigned dimLimit);
+
+  /// Adds additional local ids to the sets such that they both have the union
+  /// of the local ids in each set, without changing the set of points that
+  /// lie in `this` and `other`. The ordering of the local ids in the
+  /// sets may also be changed. After merging, if the `i^th` local variable in
+  /// one set has a known division representation, then the `i^th` local
+  /// variable in the other set either has the same division representation or
+  /// no known division representation.
+  ///
+  /// The number of dimensions and symbol ids in `this` and `other` should
+  /// match.
+  void mergeLocalIds(IntegerPolyhedron &other);
+
+  void print(raw_ostream &os) const;
+  void dump() const;
+
 protected:
+  /// Returns false if the fields corresponding to various identifier counts, or
+  /// equality/inequality buffer sizes aren't consistent; true otherwise. This
+  /// is meant to be used within an assert internally.
+  virtual bool hasConsistentState() const;
+
+  /// Checks all rows of equality/inequality constraints for trivial
+  /// contradictions (for example: 1 == 0, 0 >= 1), which may have surfaced
+  /// after elimination. Returns true if an invalid constraint is found;
+  /// false otherwise.
+  bool hasInvalidConstraint() const;
+
+  /// Returns the constant lower bound bound if isLower is true, and the upper
+  /// bound if isLower is false.
+  template <bool isLower>
+  Optional<int64_t> computeConstantLowerOrUpperBound(unsigned pos);
+
+  /// Eliminates a single identifier at `position` from equality and inequality
+  /// constraints. Returns `success` if the identifier was eliminated, and
+  /// `failure` otherwise.
+  inline LogicalResult gaussianEliminateId(unsigned position) {
+    return success(gaussianEliminateIds(position, position + 1) == 1);
+  }
+
+  /// Removes local variables using equalities. Each equality is checked if it
+  /// can be reduced to the form: `e = affine-expr`, where `e` is a local
+  /// variable and `affine-expr` is an affine expression not containing `e`.
+  /// If an equality satisfies this form, the local variable is replaced in
+  /// each constraint and then removed. The equality used to replace this local
+  /// variable is also removed.
+  void removeRedundantLocalVars();
+
+  /// Eliminates identifiers from equality and inequality constraints
+  /// in column range [posStart, posLimit).
+  /// Returns the number of variables eliminated.
+  unsigned gaussianEliminateIds(unsigned posStart, unsigned posLimit);
+
+  /// Eliminates the identifier at the specified position using Fourier-Motzkin
+  /// variable elimination, but uses Gaussian elimination if there is an
+  /// equality involving that identifier. If the result of the elimination is
+  /// integer exact, `*isResultIntegerExact` is set to true. If `darkShadow` is
+  /// set to true, a potential under approximation (subset) of the rational
+  /// shadow / exact integer shadow is computed.
+  // See implementation comments for more details.
+  virtual void fourierMotzkinEliminate(unsigned pos, bool darkShadow = false,
+                                       bool *isResultIntegerExact = nullptr);
+
+  /// Tightens inequalities given that we are dealing with integer spaces. This
+  /// is similar to the GCD test but applied to inequalities. The constant term
+  /// can be reduced to the preceding multiple of the GCD of the coefficients,
+  /// i.e.,
+  ///  64*i - 100 >= 0  =>  64*i - 128 >= 0 (since 'i' is an integer). This is a
+  /// fast method (linear in the number of coefficients).
+  void gcdTightenInequalities();
+
+  /// Normalized each constraints by the GCD of its coefficients.
+  void normalizeConstraintsByGCD();
+
+  /// Searches for a constraint with a non-zero coefficient at `colIdx` in
+  /// equality (isEq=true) or inequality (isEq=false) constraints.
+  /// Returns true and sets row found in search in `rowIdx`, false otherwise.
+  bool findConstraintWithNonZeroAt(unsigned colIdx, bool isEq,
+                                   unsigned *rowIdx) const;
+
+  /// Returns true if the pos^th column is all zero for both inequalities and
+  /// equalities.
+  bool isColZero(unsigned pos) const;
+
+  /// Prints the number of constraints, dimensions, symbols, and locals.
+  virtual void printSpace(raw_ostream &os) const;
+
+  /// A parameter that controls detection of an unrealistic number of
+  /// constraints. If the number of constraints is this many times the number of
+  /// variables, we consider such a system out of line with the intended use
+  /// case of IntegerPolyhedron.
+  // The rationale for 32 is that in the typical simplest of cases, an
+  // identifier is expected to have one lower bound and one upper bound
+  // constraint. With a level of tiling or a connection to another identifier
+  // through a div or mod, an extra pair of bounds gets added. As a limit, we
+  // don't expect an identifier to have more than 32 lower/upper/equality
+  // constraints. This is conservatively set low and can be raised if needed.
+  constexpr static unsigned kExplosionFactor = 32;
+
   /// Return the index at which the specified kind of id starts.
   unsigned getIdKindOffset(IdKind kind) const;
 
   /// Get the number of ids of the specified kind.
   unsigned getNumIdKind(IdKind kind) const;
 
   /// Removes identifiers in the column range [idStart, idLimit), and copies any
   /// remaining valid data into place, updates member variables, and resizes

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.
+// Support for linear transforms and applying them to IntegerPolyhedron.
 //
 //===----------------------------------------------------------------------===//
 
 #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 IntegerPolyhedron;
+
 class LinearTransform {
 public:
   explicit LinearTransform(Matrix &&oMatrix);
   explicit LinearTransform(const Matrix &oMatrix);
 
   // 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<unsigned, LinearTransform>
   makeTransformToColumnEchelon(Matrix m);
 
-  // Returns a FlatAffineConstraints having a constraint vector vT for every
+  // Returns a IntegerPolyhedron having a constraint vector vT for every
   // constraint vector v in fac, where T is this transform.
-  FlatAffineConstraints applyTo(const FlatAffineConstraints &fac) const;
+  IntegerPolyhedron applyTo(const IntegerPolyhedron &fac) const;
 
   // The given vector is interpreted as a row vector v. Post-multiply v with
   // this transform, say T, and return vT.
   SmallVector<int64_t, 8> postMultiplyRow(ArrayRef<int64_t> rowVec) const;
 
   // The given vector is interpreted as a column vector v. Pre-multiply v with
   // this transform, say T, and return Tv.
   SmallVector<int64_t, 8> preMultiplyColumn(ArrayRef<int64_t> colVec) const;
 
 private:
   Matrix matrix;
 };
 
 } // namespace mlir
 #endif // MLIR_ANALYSIS_LINEARTRANSFORM_H

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

 //===- Simplex.h - MLIR Simplex 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
 //
 //===----------------------------------------------------------------------===//
 //
-// Functionality to perform analysis on FlatAffineConstraints. In particular,
+// Functionality to perform analysis on IntegerPolyhedron. In particular,
 // support for performing emptiness checks and redundancy checks.
 //
 //===----------------------------------------------------------------------===//
 
 #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/IntegerPolyhedron.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;
 
 /// The Simplex class implements a version of the Simplex and Generalized Basis
 /// Reduction algorithms, which can perform analysis of integer sets with affine
 /// inequalities and equalities. A Simplex can be constructed
 /// by specifying the dimensionality of the set. It supports adding affine
 /// inequalities and equalities, and can perform emptiness checks, i.e., it can
 /// find a solution to the set of constraints if one exists, or say that the
 /// set is empty if no solution exists. Currently, this only works for bounded
 /// sets. Furthermore, it can find a subset of these constraints that are
 /// redundant, i.e. a subset of constraints that doesn't constrain the affine
 /// set further after adding the non-redundant constraints. Simplex can also be
-/// constructed from a FlatAffineConstraints object.
+/// constructed from a IntegerPolyhedron object.
 ///
 /// The implementation of the Simplex and SimplexBase classes, other than the
 /// functionality for sampling, is based on the paper
 /// "Simplify: A Theorem Prover for Program Checking"
 /// by D. Detlefs, G. Nelson, J. B. Saxe.
 ///
 /// We define variables, constraints, and unknowns. Consider the example of a
 /// two-dimensional set defined by 1 + 2x + 3y >= 0 and 2x - 3y >= 0. Here,

 /// in `Simplex`.
 
 class SimplexBase {
 public:
   enum class Direction { Up, Down };
 
   SimplexBase() = delete;
   explicit SimplexBase(unsigned nVar);
-  explicit SimplexBase(const FlatAffineConstraints &constraints);
+  explicit SimplexBase(const IntegerPolyhedron &constraints);
 
   /// Returns true if the tableau is empty (has conflicting constraints),
   /// false otherwise.
   bool isEmpty() const;
 
   /// Add an inequality to the tableau. If coeffs is c_0, c_1, ... c_n, where n
   /// is the current number of variables, then the corresponding inequality is
   /// c_n + c_0*x_0 + c_1*x_1 + ... + c_{n-1}*x_{n-1} >= 0.
public:
   void markEmpty();
 
   /// Get a snapshot of the current state. This is used for rolling back.
   unsigned getSnapshot() const;
 
   /// 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);
+  /// Add all the constraints from the given IntegerPolyhedron.
+  void intersectIntegerPolyhedron(const IntegerPolyhedron &fac);
 
   /// Returns a rational sample point. This should not be called when Simplex is
   /// empty.
   SmallVector<Fraction, 8> getRationalSample() const;
 
   /// Returns the current sample point if it is integral. Otherwise, returns
   /// None.
   Optional<SmallVector<int64_t, 8>> getSamplePointIfIntegral() const;
protected:
   /// These hold information about each unknown.
   SmallVector<Unknown, 8> con, var;
 };
 
 class Simplex : public SimplexBase {
 public:
   Simplex() = delete;
   explicit Simplex(unsigned nVar) : SimplexBase(nVar) {}
-  explicit Simplex(const FlatAffineConstraints &constraints)
+  explicit Simplex(const IntegerPolyhedron &constraints)
       : SimplexBase(constraints) {}
 
   /// Compute the maximum or minimum value of the given row, depending on
   /// direction. The specified row is never pivoted. On return, the row may
   /// have a negative sample value if the direction is down.
   ///
   /// Returns a Fraction denoting the optimum, or a null value if no optimum
   /// exists, i.e., if the expression is unbounded in this direction.
public:
   /// Check if the specified inequality already holds in the polytope.
   bool isRedundantInequality(ArrayRef<int64_t> coeffs);
 
   /// Check if the specified equality already holds in the polytope.
   bool isRedundantEquality(ArrayRef<int64_t> coeffs);
 
   /// Returns true if this Simplex's polytope is a rational subset of `fac`.
   /// Otherwise, returns false.
-  bool isRationalSubsetOf(const FlatAffineConstraints &fac);
+  bool isRationalSubsetOf(const IntegerPolyhedron &fac);
 
 private:
   friend class GBRSimplex;
 
   /// 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.

mlir/include/mlir/Analysis/PresburgerSet.h

 //===- Set.h - MLIR PresburgerSet 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
 //
 //===----------------------------------------------------------------------===//
 //
-// A class to represent unions of FlatAffineConstraints.
+// A class to represent unions of IntegerPolyhedron.
 //
 //===----------------------------------------------------------------------===//
 
 #ifndef MLIR_ANALYSIS_PRESBURGERSET_H
 #define MLIR_ANALYSIS_PRESBURGERSET_H
 
-#include "mlir/Analysis/AffineStructures.h"
+#include "mlir/Analysis/Presburger/IntegerPolyhedron.h"
 
 namespace mlir {
 
-/// This class can represent a union of FlatAffineConstraints, with support for
+/// This class can represent a union of IntegerPolyhedron, with support for
 /// union, intersection, subtraction and complement operations, as well as
 /// sampling.
 ///
-/// The FlatAffineConstraints (FACs) are stored in a vector, and the set
+/// The IntegerPolyhedron (FACs) are stored in a vector, and the set
 /// represents the union of these FACs. An empty list corresponds to the empty
 /// set.
 ///
 /// Note that there are no invariants guaranteed on the list of FACs other than
 /// that they are all in the same space, i.e., they all have the same number of
 /// dimensions and symbols. For example, the FACs may overlap each other.
 class PresburgerSet {
 public:
-  explicit PresburgerSet(const FlatAffineConstraints &fac);
+  explicit PresburgerSet(const IntegerPolyhedron &fac);
 
   /// Return the number of FACs in the union.
   unsigned getNumFACs() const;
 
   /// Return the number of real dimensions.
   unsigned getNumDims() const;
 
   /// Return the number of symbolic dimensions.
   unsigned getNumSyms() const;
 
-  /// Return a reference to the list of FlatAffineConstraints.
-  ArrayRef<FlatAffineConstraints> getAllFlatAffineConstraints() const;
+  /// Return a reference to the list of IntegerPolyhedron.
+  ArrayRef<IntegerPolyhedron> getAllIntegerPolyhedron() const;
 
-  /// Return the FlatAffineConstraints at the specified index.
-  const FlatAffineConstraints &getFlatAffineConstraints(unsigned index) const;
+  /// Return the IntegerPolyhedron at the specified index.
+  const IntegerPolyhedron &getIntegerPolyhedron(unsigned index) const;
 
   /// Mutate this set, turning it into the union of this set and the given
-  /// FlatAffineConstraints.
-  void unionFACInPlace(const FlatAffineConstraints &fac);
+  /// IntegerPolyhedron.
+  void unionFACInPlace(const IntegerPolyhedron &fac);
 
   /// Mutate this set, turning it into the union of this set and the given set.
   void unionSetInPlace(const PresburgerSet &set);
 
   /// Return the union of this set and the given set.
   PresburgerSet unionSet(const PresburgerSet &set) const;
 
   /// Return the intersection of this set and the given set.
public:
   PresburgerSet coalesce() const;
 
 private:
   /// Construct an empty PresburgerSet.
   PresburgerSet(unsigned nDim = 0, unsigned nSym = 0)
       : nDim(nDim), nSym(nSym) {}
 
   /// Return the set difference fac \ set.
-  static PresburgerSet getSetDifference(FlatAffineConstraints fac,
+  static PresburgerSet getSetDifference(IntegerPolyhedron fac,
                                         const PresburgerSet &set);
 
   /// Number of identifiers corresponding to real dimensions.
   unsigned nDim;
 
   /// Number of symbolic dimensions, unknown but constant for analysis, as in
-  /// FlatAffineConstraints.
+  /// IntegerPolyhedron.
   unsigned nSym;
 
   /// The list of flatAffineConstraints that this set is the union of.
-  SmallVector<FlatAffineConstraints, 2> flatAffineConstraints;
+  SmallVector<IntegerPolyhedron, 2> flatAffineConstraints;
 };
 
 } // namespace mlir
 
 #endif // MLIR_ANALYSIS_PRESBURGERSET_H

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/LinearTransform.h"
-#include "mlir/Analysis/Presburger/Simplex.h"
 #include "mlir/Analysis/Presburger/Utils.h"
 #include "mlir/Dialect/Affine/IR/AffineOps.h"
 #include "mlir/Dialect/Affine/IR/AffineValueMap.h"
 #include "mlir/Dialect/Arithmetic/IR/Arithmetic.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/STLExtras.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;
 
 namespace {
 
 // See comments for SimpleAffineExprFlattener.
 // An AffineExprFlattener extends a SimpleAffineExprFlattener by recording
 // constraint information associated with mod's, floordiv's, and ceildiv's
 // in FlatAffineConstraints 'localVarCst'.
 struct AffineExprFlattener : public SimpleAffineExprFlattener {
void FlatAffineValueConstraints::addAffineIfOpDomain(AffineIfOp ifOp) {
 
   // Merge the constraints from ifOp to the current domain. We need first merge
   // and align the IDs from both constraints, and then append the constraints
   // from the ifOp into the current one.
   mergeAndAlignIdsWithOther(0, &cst);
   append(cst);
 }
 
-// Searches for a constraint with a non-zero coefficient at `colIdx` in
-// equality (isEq=true) or inequality (isEq=false) constraints.
-// Returns true and sets row found in search in `rowIdx`, false otherwise.
-bool FlatAffineConstraints::findConstraintWithNonZeroAt(
-    unsigned colIdx, bool isEq, unsigned *rowIdx) const {
-  assert(colIdx < getNumCols() && "position out of bounds");
-  auto at = [&](unsigned rowIdx) -> int64_t {
-    return isEq ? atEq(rowIdx, colIdx) : atIneq(rowIdx, colIdx);
-  };
-  unsigned e = isEq ? getNumEqualities() : getNumInequalities();
-  for (*rowIdx = 0; *rowIdx < e; ++(*rowIdx)) {
-    if (at(*rowIdx) != 0) {
-      return true;
-    }
-  }
-  return false;
-}
-
-// Normalizes the coefficient values across all columns in `rowIdx` by their
-// GCD in equality or inequality constraints as specified by `isEq`.
-template <bool isEq>
-static void normalizeConstraintByGCD(FlatAffineConstraints *constraints,
-                                     unsigned rowIdx) {
-  auto at = [&](unsigned colIdx) -> int64_t {
-    return isEq ? constraints->atEq(rowIdx, colIdx)
-                : constraints->atIneq(rowIdx, colIdx);
-  };
-  uint64_t gcd = std::abs(at(0));
-  for (unsigned j = 1, e = constraints->getNumCols(); j < e; ++j) {
-    gcd = llvm::GreatestCommonDivisor64(gcd, std::abs(at(j)));
-  }
-  if (gcd > 0 && gcd != 1) {
-    for (unsigned j = 0, e = constraints->getNumCols(); j < e; ++j) {
-      int64_t v = at(j) / static_cast<int64_t>(gcd);
-      isEq ? constraints->atEq(rowIdx, j) = v
-           : constraints->atIneq(rowIdx, j) = v;
-    }
-  }
-}
-
-void FlatAffineConstraints::normalizeConstraintsByGCD() {
-  for (unsigned i = 0, e = getNumEqualities(); i < e; ++i) {
-    normalizeConstraintByGCD</*isEq=*/true>(this, i);
-  }
-  for (unsigned i = 0, e = getNumInequalities(); i < e; ++i) {
-    normalizeConstraintByGCD</*isEq=*/false>(this, i);
-  }
-}
-
-bool FlatAffineConstraints::hasConsistentState() const {
-  if (!inequalities.hasConsistentState())
-    return false;
-  if (!equalities.hasConsistentState())
-    return false;
-
-  // Catches errors where numDims, numSymbols, numIds aren't consistent.
-  if (numDims > numIds || numSymbols > numIds || numDims + numSymbols > numIds)
-    return false;
-
-  return true;
-}
-
 bool FlatAffineValueConstraints::hasConsistentState() const {
   return FlatAffineConstraints::hasConsistentState() &&
          values.size() == getNumIds();
 }
 
-bool FlatAffineConstraints::hasInvalidConstraint() const {
-  assert(hasConsistentState());
-  auto check = [&](bool isEq) -> bool {
-    unsigned numCols = getNumCols();
-    unsigned numRows = isEq ? getNumEqualities() : getNumInequalities();
-    for (unsigned i = 0, e = numRows; i < e; ++i) {
-      unsigned j;
-      for (j = 0; j < numCols - 1; ++j) {
-        int64_t v = isEq ? atEq(i, j) : atIneq(i, j);
-        // Skip rows with non-zero variable coefficients.
-        if (v != 0)
-          break;
-      }
-      if (j < numCols - 1) {
-        continue;
-      }
-      // Check validity of constant term at 'numCols - 1' w.r.t 'isEq'.
-      // Example invalid constraints include: '1 == 0' or '-1 >= 0'
-      int64_t v = isEq ? atEq(i, numCols - 1) : atIneq(i, numCols - 1);
-      if ((isEq && v != 0) || (!isEq && v < 0)) {
-        return true;
-      }
-    }
-    return false;
-  };
-  if (check(/*isEq=*/true))
-    return true;
-  return check(/*isEq=*/false);
-}
-
-/// Eliminate identifier from constraint at `rowIdx` based on coefficient at
-/// pivotRow, pivotCol. Columns in range [elimColStart, pivotCol) will not be
-/// updated as they have already been eliminated.
-static void eliminateFromConstraint(FlatAffineConstraints *constraints,
-                                    unsigned rowIdx, unsigned pivotRow,
-                                    unsigned pivotCol, unsigned elimColStart,
-                                    bool isEq) {
-  // Skip if equality 'rowIdx' if same as 'pivotRow'.
-  if (isEq && rowIdx == pivotRow)
-    return;
-  auto at = [&](unsigned i, unsigned j) -> int64_t {
-    return isEq ? constraints->atEq(i, j) : constraints->atIneq(i, j);
-  };
-  int64_t leadCoeff = at(rowIdx, pivotCol);
-  // Skip if leading coefficient at 'rowIdx' is already zero.
-  if (leadCoeff == 0)
-    return;
-  int64_t pivotCoeff = constraints->atEq(pivotRow, pivotCol);
-  int64_t sign = (leadCoeff * pivotCoeff > 0) ? -1 : 1;
-  int64_t lcm = mlir::lcm(pivotCoeff, leadCoeff);
-  int64_t pivotMultiplier = sign * (lcm / std::abs(pivotCoeff));
-  int64_t rowMultiplier = lcm / std::abs(leadCoeff);
-
-  unsigned numCols = constraints->getNumCols();
-  for (unsigned j = 0; j < numCols; ++j) {
-    // Skip updating column 'j' if it was just eliminated.
-    if (j >= elimColStart && j < pivotCol)
-      continue;
-    int64_t v = pivotMultiplier * constraints->atEq(pivotRow, j) +
-                rowMultiplier * at(rowIdx, j);
-    isEq ? constraints->atEq(rowIdx, j) = v
-         : constraints->atIneq(rowIdx, j) = v;
-  }
-}
-
 void FlatAffineValueConstraints::removeIdRange(unsigned idStart,
                                                unsigned idLimit) {
   FlatAffineConstraints::removeIdRange(idStart, idLimit);
   values.erase(values.begin() + idStart, values.begin() + idLimit);
 }
 
-/// Returns the position of the identifier that has the minimum <number of lower
-/// bounds> times <number of upper bounds> from the specified range of
-/// identifiers [start, end). It is often best to eliminate in the increasing
-/// order of these counts when doing Fourier-Motzkin elimination since FM adds
-/// that many new constraints.
-static unsigned getBestIdToEliminate(const FlatAffineConstraints &cst,
-                                     unsigned start, unsigned end) {
-  assert(start < cst.getNumIds() && end < cst.getNumIds() + 1);
-
-  auto getProductOfNumLowerUpperBounds = [&](unsigned pos) {
-    unsigned numLb = 0;
-    unsigned numUb = 0;
-    for (unsigned r = 0, e = cst.getNumInequalities(); r < e; r++) {
-      if (cst.atIneq(r, pos) > 0) {
-        ++numLb;
-      } else if (cst.atIneq(r, pos) < 0) {
-        ++numUb;
-      }
-    }
-    return numLb * numUb;
-  };
-
-  unsigned minLoc = start;
-  unsigned min = getProductOfNumLowerUpperBounds(start);
-  for (unsigned c = start + 1; c < end; c++) {
-    unsigned numLbUbProduct = getProductOfNumLowerUpperBounds(c);
-    if (numLbUbProduct < min) {
-      min = numLbUbProduct;
-      minLoc = c;
-    }
-  }
-  return minLoc;
-}
-
-// Checks for emptiness of the set by eliminating identifiers successively and
-// using the GCD test (on all equality constraints) and checking for trivially
-// invalid constraints. Returns 'true' if the constraint system is found to be
-// empty; false otherwise.
-bool FlatAffineConstraints::isEmpty() const {
-  if (isEmptyByGCDTest() || hasInvalidConstraint())
-    return true;
-
-  FlatAffineConstraints tmpCst(*this);
-
-  // First, eliminate as many local variables as possible using equalities.
-  tmpCst.removeRedundantLocalVars();
-  if (tmpCst.isEmptyByGCDTest() || tmpCst.hasInvalidConstraint())
-    return true;
-
-  // Eliminate as many identifiers as possible using Gaussian elimination.
-  unsigned currentPos = 0;
-  while (currentPos < tmpCst.getNumIds()) {
-    tmpCst.gaussianEliminateIds(currentPos, tmpCst.getNumIds());
-    ++currentPos;
-    // We check emptiness through trivial checks after eliminating each ID to
-    // detect emptiness early. Since the checks isEmptyByGCDTest() and
-    // hasInvalidConstraint() are linear time and single sweep on the constraint
-    // buffer, this appears reasonable - but can optimize in the future.
-    if (tmpCst.hasInvalidConstraint() || tmpCst.isEmptyByGCDTest())
-      return true;
-  }
-
-  // Eliminate the remaining using FM.
-  for (unsigned i = 0, e = tmpCst.getNumIds(); i < e; i++) {
-    tmpCst.fourierMotzkinEliminate(
-        getBestIdToEliminate(tmpCst, 0, tmpCst.getNumIds()));
-    // Check for a constraint explosion. This rarely happens in practice, but
-    // this check exists as a safeguard against improperly constructed
-    // constraint systems or artificially created arbitrarily complex systems
-    // that aren't the intended use case for FlatAffineConstraints. This is
-    // needed since FM has a worst case exponential complexity in theory.
-    if (tmpCst.getNumConstraints() >= kExplosionFactor * getNumIds()) {
-      LLVM_DEBUG(llvm::dbgs() << "FM constraint explosion detected\n");
-      return false;
-    }
-
-    // FM wouldn't have modified the equalities in any way. So no need to again
-    // run GCD test. Check for trivial invalid constraints.
-    if (tmpCst.hasInvalidConstraint())
-      return true;
-  }
-  return false;
-}
-
-// Runs the GCD test on all equality constraints. Returns 'true' if this test
-// fails on any equality. Returns 'false' otherwise.
-// This test can be used to disprove the existence of a solution. If it returns
-// true, no integer solution to the equality constraints can exist.
-//
-// GCD test definition:
-//
-// The equality constraint:
-//
-//  c_1*x_1 + c_2*x_2 + ... + c_n*x_n = c_0
-//
-// has an integer solution iff:
-//
-//  GCD of c_1, c_2, ..., c_n divides c_0.
-//
-bool FlatAffineConstraints::isEmptyByGCDTest() const {
-  assert(hasConsistentState());
-  unsigned numCols = getNumCols();
-  for (unsigned i = 0, e = getNumEqualities(); i < e; ++i) {
-    uint64_t gcd = std::abs(atEq(i, 0));
-    for (unsigned j = 1; j < numCols - 1; ++j) {
-      gcd = llvm::GreatestCommonDivisor64(gcd, std::abs(atEq(i, j)));
-    }
-    int64_t v = std::abs(atEq(i, numCols - 1));
-    if (gcd > 0 && (v % gcd != 0)) {
-      return true;
-    }
-  }
-  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.
-//
-// 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.
-  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.getNumIds(), 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.getNumIds(), 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);
-}
-
-bool FlatAffineConstraints::isIntegerEmpty() const {
-  return !findIntegerSample().hasValue();
-}
-
-/// 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 infinity 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 B by removing all constraints that involve
-/// the unbounded dimensions and then deleting the unbounded dimensions. Note
-/// that B is a Bounded set.
-///
-/// 3) Try to obtain a sample from B using the GBR sampling
-/// algorithm. If no sample is found, return that S is empty.
-///
-/// 4) Otherwise, substitute the obtained sample into S*T to obtain a set
-/// C. C is a full-dimensional Cone and always contains a sample.
-///
-/// 5) Obtain an integer sample from C.
-///
-/// 6) Return T*v, where v is the concatenation of the samples from B and C.
-///
-/// The following is a sketch of a proof that
-/// a) If the algorithm returns empty, then S is empty.
-/// b) If the algorithm returns a sample, it is a valid sample in S.
-///
-/// The algorithm returns empty only if B is empty, in which case S*T is
-/// certainly empty since B was obtained by removing constraints and then
-/// deleting unconstrained dimensions from S*T. Since T is unimodular, a vector
-/// v is in S*T iff T*v is in S. So in this case, since
-/// S*T is empty, S is empty too.
-///
-/// Otherwise, the algorithm substitutes the sample from B into S*T. All the
-/// constraints of S*T that did not involve unbounded dimensions are satisfied
-/// by this substitution. All dimensions in the linear span of the dimensions
-/// outside the prefix are unbounded in S*T (step 1). Substituting values for
-/// the bounded dimensions cannot make these dimensions bounded, and these are
-/// the only remaining dimensions in C, so C is unbounded along every vector (in
-/// the positive or negative direction, or both). C is hence a full-dimensional
-/// cone and therefore always contains an integer point.
-///
-/// Concatenating the samples from B and C gives a sample v in S*T, so the
-/// returned sample T*v is a sample in S.
-Optional<SmallVector<int64_t, 8>>
-FlatAffineConstraints::findIntegerSample() const {
-  // First, try the GCD test heuristic.
-  if (isEmptyByGCDTest())
-    return {};
-
-  Simplex simplex(*this);
-  if (simplex.isEmpty())
-    return {};
-
-  // For a bounded set, we directly call into the GBR sampling algorithm.
-  if (!simplex.isUnbounded())
-    return simplex.findIntegerSample();
-
-  // 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<unsigned, LinearTransform> result =
-      LinearTransform::makeTransformToColumnEchelon(std::move(m));
-  const LinearTransform &transform = result.second;
-  // 1) Apply T to S to obtain S*T.
-  FlatAffineConstraints transformedSet = transform.applyTo(*this);
-
-  // 2) Remove the unbounded dimensions and constraints involving them to
-  // obtain a bounded set.
-  FlatAffineConstraints boundedSet = transformedSet;
-  unsigned numBoundedDims = result.first;
-  unsigned numUnboundedDims = getNumIds() - numBoundedDims;
-  removeConstraintsInvolvingSuffixDims(boundedSet, numUnboundedDims);
-  boundedSet.removeIdRange(numBoundedDims, boundedSet.getNumIds());
-
-  // 3) Try to obtain a sample from the bounded set.
-  Optional<SmallVector<int64_t, 8>> boundedSample =
-      Simplex(boundedSet).findIntegerSample();
-  if (!boundedSample)
-    return {};
-  assert(boundedSet.containsPoint(*boundedSample) &&
-         "Simplex returned an invalid sample!");
-
-  // 4) Substitute the values of the bounded dimensions into S*T to obtain a
-  // full-dimensional cone, which necessarily contains an integer sample.
-  transformedSet.setAndEliminate(0, *boundedSample);
-  FlatAffineConstraints &cone = transformedSet;
-
-  // 5) Obtain an integer sample from the cone.
-  //
-  // We shrink the cone such that for any rational point in the shrunken cone,
-  // rounding up each of the point's coordinates produces a point that still
-  // lies in the original cone.
-  //
-  // Rounding up a point x adds a number e_i in [0, 1) to each coordinate x_i.
-  // For each inequality sum_i a_i x_i + c >= 0 in the original cone, the
-  // shrunken cone will have the inequality tightened by some amount s, such
-  // that if x satisfies the shrunken cone's tightened inequality, then x + e
-  // satisfies the original inequality, i.e.,
-  //
-  // sum_i a_i x_i + c + s >= 0 implies sum_i a_i (x_i + e_i) + c >= 0
-  //
-  // for any e_i values in [0, 1). In fact, we will handle the slightly more
-  // general case where e_i can be in [0, 1]. For example, consider the
-  // inequality 2x_1 - 3x_2 - 7x_3 - 6 >= 0, and let x = (3, 0, 0). How low
-  // could the LHS go if we added a number in [0, 1] to each coordinate? The LHS
-  // is minimized when we add 1 to the x_i with negative coefficient a_i and
-  // keep the other x_i the same. In the example, we would get x = (3, 1, 1),
-  // changing the value of the LHS by -3 + -7 = -10.
-  //
-  // In general, the value of the LHS can change by at most the sum of the
-  // negative a_i, so we accomodate this by shifting the inequality by this
-  // amount for the shrunken cone.
-  for (unsigned i = 0, e = cone.getNumInequalities(); i < e; ++i) {
-    for (unsigned j = 0; j < cone.numIds; ++j) {
-      int64_t coeff = cone.atIneq(i, j);
-      if (coeff < 0)
-        cone.atIneq(i, cone.numIds) += coeff;
-    }
-  }
-
-  // Obtain an integer sample in the cone by rounding up a rational point from
-  // the shrunken cone. Shrinking the cone amounts to shifting its apex
-  // "inwards" without changing its "shape"; the shrunken cone is still a
-  // full-dimensional cone and is hence non-empty.
-  Simplex shrunkenConeSimplex(cone);
-  assert(!shrunkenConeSimplex.isEmpty() && "Shrunken cone cannot be empty!");
-  SmallVector<Fraction, 8> shrunkenConeSample =
-      shrunkenConeSimplex.getRationalSample();
-
-  SmallVector<int64_t, 8> coneSample(llvm::map_range(shrunkenConeSample, ceil));
-
-  // 6) Return transform * concat(boundedSample, coneSample).
-  SmallVector<int64_t, 8> &sample = boundedSample.getValue();
-  sample.append(coneSample.begin(), coneSample.end());
-  return transform.preMultiplyColumn(sample);
-}
-
-/// Helper to evaluate an affine expression at a point.
-/// The expression is a list of coefficients for the dimensions followed by the
-/// constant term.
-static int64_t valueAt(ArrayRef<int64_t> expr, ArrayRef<int64_t> point) {
-  assert(expr.size() == 1 + point.size() &&
-         "Dimensionalities of point and expression don't match!");
-  int64_t value = expr.back();
-  for (unsigned i = 0; i < point.size(); ++i)
-    value += expr[i] * point[i];
-  return value;
-}
-
-/// A point satisfies an equality iff the value of the equality at the
-/// expression is zero, and it satisfies an inequality iff the value of the
-/// inequality at that point is non-negative.
-bool FlatAffineConstraints::containsPoint(ArrayRef<int64_t> point) const {
-  for (unsigned i = 0, e = getNumEqualities(); i < e; ++i) {
-    if (valueAt(getEquality(i), point) != 0)
-      return false;
-  }
-  for (unsigned i = 0, e = getNumInequalities(); i < e; ++i) {
-    if (valueAt(getInequality(i), point) < 0)
-      return false;
-  }
-  return true;
-}
-
-void FlatAffineConstraints::getLocalReprs(
-    std::vector<llvm::Optional<std::pair<unsigned, unsigned>>> &repr) const {
-  std::vector<SmallVector<int64_t, 8>> dividends(getNumLocalIds());
-  SmallVector<unsigned, 4> denominators(getNumLocalIds());
-  getLocalReprs(dividends, denominators, repr);
-}
-
-void FlatAffineConstraints::getLocalReprs(
-    std::vector<SmallVector<int64_t, 8>> &dividends,
-    SmallVector<unsigned, 4> &denominators) const {
-  std::vector<llvm::Optional<std::pair<unsigned, unsigned>>> repr(
-      getNumLocalIds());
-  getLocalReprs(dividends, denominators, repr);
-}
-
-void FlatAffineConstraints::getLocalReprs(
-    std::vector<SmallVector<int64_t, 8>> &dividends,
-    SmallVector<unsigned, 4> &denominators,
-    std::vector<llvm::Optional<std::pair<unsigned, unsigned>>> &repr) const {
-
-  repr.resize(getNumLocalIds());
-  dividends.resize(getNumLocalIds());
-  denominators.resize(getNumLocalIds());
-
-  SmallVector<bool, 8> foundRepr(getNumIds(), false);
-  for (unsigned i = 0, e = getNumDimAndSymbolIds(); i < e; ++i)
-    foundRepr[i] = true;
-
-  unsigned divOffset = getNumDimAndSymbolIds();
-  bool changed;
-  do {
-    // Each time changed is true, at end of this iteration, one or more local
-    // vars have been detected as floor divs.
-    changed = false;
-    for (unsigned i = 0, e = getNumLocalIds(); i < e; ++i) {
-      if (!foundRepr[i + divOffset]) {
-        if (auto res = presburger_utils::computeSingleVarRepr(
-                *this, foundRepr, divOffset + i, dividends[i],
-                denominators[i])) {
-          foundRepr[i + divOffset] = true;
-          repr[i] = res;
-          changed = true;
-        }
-      }
-    }
-  } while (changed);
-
-  // Set 0 denominator for identifiers for which no division representation
-  // could be found.
-  for (unsigned i = 0, e = repr.size(); i < e; ++i)
-    if (!repr[i].hasValue())
-      denominators[i] = 0;
-}
-
-/// Tightens inequalities given that we are dealing with integer spaces. This is
-/// analogous to the GCD test but applied to inequalities. The constant term can
-/// be reduced to the preceding multiple of the GCD of the coefficients, i.e.,
-///  64*i - 100 >= 0  =>  64*i - 128 >= 0 (since 'i' is an integer). This is a
-/// fast method - linear in the number of coefficients.
-// Example on how this affects practical cases: consider the scenario:
-// 64*i >= 100, j = 64*i; without a tightening, elimination of i would yield
-// j >= 100 instead of the tighter (exact) j >= 128.
-void FlatAffineConstraints::gcdTightenInequalities() {
-  unsigned numCols = getNumCols();
-  for (unsigned i = 0, e = getNumInequalities(); i < e; ++i) {
-    uint64_t gcd = std::abs(atIneq(i, 0));
-    for (unsigned j = 1; j < numCols - 1; ++j) {
-      gcd = llvm::GreatestCommonDivisor64(gcd, std::abs(atIneq(i, j)));
-    }
-    if (gcd > 0 && gcd != 1) {
-      int64_t gcdI = static_cast<int64_t>(gcd);
-      // Tighten the constant term and normalize the constraint by the GCD.
-      atIneq(i, numCols - 1) = mlir::floorDiv(atIneq(i, numCols - 1), gcdI);
-      for (unsigned j = 0, e = numCols - 1; j < e; ++j)
-        atIneq(i, j) /= gcdI;
-    }
-  }
-}
-
-// Eliminates all identifier variables in column range [posStart, posLimit).
-// Returns the number of variables eliminated.
-unsigned FlatAffineConstraints::gaussianEliminateIds(unsigned posStart,
-                                                     unsigned posLimit) {
-  // Return if identifier positions to eliminate are out of range.
-  assert(posLimit <= numIds);
-  assert(hasConsistentState());
-
-  if (posStart >= posLimit)
-    return 0;
-
-  gcdTightenInequalities();
-
-  unsigned pivotCol = 0;
-  for (pivotCol = posStart; pivotCol < posLimit; ++pivotCol) {
-    // Find a row which has a non-zero coefficient in column 'j'.
-    unsigned pivotRow;
-    if (!findConstraintWithNonZeroAt(pivotCol, /*isEq=*/true, &pivotRow)) {
-      // No pivot row in equalities with non-zero at 'pivotCol'.
-      if (!findConstraintWithNonZeroAt(pivotCol, /*isEq=*/false, &pivotRow)) {
-        // If inequalities are also non-zero in 'pivotCol', it can be
-        // eliminated.
-        continue;
-      }
-      break;
-    }
-
-    // Eliminate identifier at 'pivotCol' from each equality row.
-    for (unsigned i = 0, e = getNumEqualities(); i < e; ++i) {
-      eliminateFromConstraint(this, i, pivotRow, pivotCol, posStart,
-                              /*isEq=*/true);
-      normalizeConstraintByGCD</*isEq=*/true>(this, i);
-    }
-
-    // Eliminate identifier at 'pivotCol' from each inequality row.
-    for (unsigned i = 0, e = getNumInequalities(); i < e; ++i) {
-      eliminateFromConstraint(this, i, pivotRow, pivotCol, posStart,
-                              /*isEq=*/false);
-      normalizeConstraintByGCD</*isEq=*/false>(this, i);
-    }
-    removeEquality(pivotRow);
-    gcdTightenInequalities();
-  }
-  // Update position limit based on number eliminated.
-  posLimit = pivotCol;
-  // Remove eliminated columns from all constraints.
-  removeIdRange(posStart, posLimit);
-  return posLimit - posStart;
-}
-
 // Determine whether the identifier at 'pos' (say id_r) can be expressed as
 // modulo of another known identifier (say id_n) w.r.t a constant. For example,
 // if the following constraints hold true:
 // ```
 // 0 <= id_r <= divisor - 1
 // id_n - (divisor * q_expr) = id_r
 // ```
 // where `id_n` is a known identifier (called dividend), and `q_expr` is an
for (unsigned c = 0, f = cst.getNumIds(); c < f; c++)
     if (dividend[c] != 0)
       dividendExpr = dividendExpr + dividend[c] * exprs[c];
 
   // Successfully detected the floordiv.
   exprs[pos] = dividendExpr.floorDiv(divisor);
   return true;
 }
 
-// Fills an inequality row with the value 'val'.
-static inline void fillInequality(FlatAffineConstraints *cst, unsigned r,
-                                  int64_t val) {
-  for (unsigned c = 0, f = cst->getNumCols(); c < f; c++) {
-    cst->atIneq(r, c) = val;
-  }
-}
-
-// Negates an inequality.
-static inline void negateInequality(FlatAffineConstraints *cst, unsigned r) {
-  for (unsigned c = 0, f = cst->getNumCols(); c < f; c++) {
-    cst->atIneq(r, c) = -cst->atIneq(r, c);
-  }
-}
-
-// A more complex check to eliminate redundant inequalities. Uses FourierMotzkin
-// to check if a constraint is redundant.
-void FlatAffineConstraints::removeRedundantInequalities() {
-  SmallVector<bool, 32> redun(getNumInequalities(), false);
-  // To check if an inequality is redundant, we replace the inequality by its
-  // complement (for eg., i - 1 >= 0 by i <= 0), and check if the resulting
-  // system is empty. If it is, the inequality is redundant.
-  FlatAffineConstraints tmpCst(*this);
-  for (unsigned r = 0, e = getNumInequalities(); r < e; r++) {
-    // Change the inequality to its complement.
-    negateInequality(&tmpCst, r);
-    tmpCst.atIneq(r, tmpCst.getNumCols() - 1)--;
-    if (tmpCst.isEmpty()) {
-      redun[r] = true;
-      // Zero fill the redundant inequality.
-      fillInequality(this, r, /*val=*/0);
-      fillInequality(&tmpCst, r, /*val=*/0);
-    } else {
-      // Reverse the change (to avoid recreating tmpCst each time).
-      tmpCst.atIneq(r, tmpCst.getNumCols() - 1)++;
-      negateInequality(&tmpCst, r);
-    }
-  }
-
-  // Scan to get rid of all rows marked redundant, in-place.
-  auto copyRow = [&](unsigned src, unsigned dest) {
-    if (src == dest)
-      return;
-    for (unsigned c = 0, e = getNumCols(); c < e; c++) {
-      atIneq(dest, c) = atIneq(src, c);
-    }
-  };
-  unsigned pos = 0;
-  for (unsigned r = 0, e = getNumInequalities(); r < e; r++) {
-    if (!redun[r])
-      copyRow(r, pos++);
-  }
-  inequalities.resizeVertically(pos);
-}
-
-// A more complex check to eliminate redundant inequalities and equalities. Uses
-// Simplex to check if a constraint is redundant.
-void FlatAffineConstraints::removeRedundantConstraints() {
-  // First, we run gcdTightenInequalities. This allows us to catch some
-  // constraints which are not redundant when considering rational solutions
-  // but are redundant in terms of integer solutions.
-  gcdTightenInequalities();
-  Simplex simplex(*this);
-  simplex.detectRedundant();
-
-  auto copyInequality = [&](unsigned src, unsigned dest) {
-    if (src == dest)
-      return;
-    for (unsigned c = 0, e = getNumCols(); c < e; c++)
-      atIneq(dest, c) = atIneq(src, c);
-  };
-  unsigned pos = 0;
-  unsigned numIneqs = getNumInequalities();
-  // Scan to get rid of all inequalities marked redundant, in-place. In Simplex,
-  // the first constraints added are the inequalities.
-  for (unsigned r = 0; r < numIneqs; r++) {
-    if (!simplex.isMarkedRedundant(r))
-      copyInequality(r, pos++);
-  }
-  inequalities.resizeVertically(pos);
-
-  // Scan to get rid of all equalities marked redundant, in-place. In Simplex,
-  // after the inequalities, a pair of constraints for each equality is added.
-  // An equality is redundant if both the inequalities in its pair are
-  // redundant.
-  auto copyEquality = [&](unsigned src, unsigned dest) {
-    if (src == dest)
-      return;
-    for (unsigned c = 0, e = getNumCols(); c < e; c++)
-      atEq(dest, c) = atEq(src, c);
-  };
-  pos = 0;
-  for (unsigned r = 0, e = getNumEqualities(); r < e; r++) {
-    if (!(simplex.isMarkedRedundant(numIneqs + 2 * r) &&
-          simplex.isMarkedRedundant(numIneqs + 2 * r + 1)))
-      copyEquality(r, pos++);
-  }
-  equalities.resizeVertically(pos);
-}
-
-/// Eliminate `pos2^th` local identifier, replacing its every instance with
-/// `pos1^th` local identifier. This function is intended to be used to remove
-/// redundancy when local variables at position `pos1` and `pos2` are restricted
-/// to have the same value.
-static void eliminateRedundantLocalId(FlatAffineConstraints &fac, unsigned pos1,
-                                      unsigned pos2) {
-
-  assert(pos1 < fac.getNumLocalIds() && "Invalid local id position");
-  assert(pos2 < fac.getNumLocalIds() && "Invalid local id position");
-
-  unsigned localOffset = fac.getNumDimAndSymbolIds();
-  pos1 += localOffset;
-  pos2 += localOffset;
-  for (unsigned i = 0, e = fac.getNumInequalities(); i < e; ++i)
-    fac.atIneq(i, pos1) += fac.atIneq(i, pos2);
-  for (unsigned i = 0, e = fac.getNumEqualities(); i < e; ++i)
-    fac.atEq(i, pos1) += fac.atEq(i, pos2);
-  fac.removeId(pos2);
-}
-
-/// Adds additional local ids to the sets such that they both have the union
-/// of the local ids in each set, without changing the set of points that
-/// lie in `this` and `other`.
-///
-/// To detect local ids that always take the same in both sets, each local id is
-/// represented as a floordiv with constant denominator in terms of other ids.
-/// After extracting these divisions, local ids with the same division
-/// representation are considered duplicate and are merged. It is possible that
-/// division representation for some local id cannot be obtained, and thus these
-/// local ids are not considered for detecting duplicates.
-void FlatAffineConstraints::mergeLocalIds(FlatAffineConstraints &other) {
-  assert(getNumDimIds() == other.getNumDimIds() &&
-         "Number of dimension ids should match");
-  assert(getNumSymbolIds() == other.getNumSymbolIds() &&
-         "Number of symbol ids should match");
-
-  FlatAffineConstraints &facA = *this;
-  FlatAffineConstraints &facB = other;
-
-  // Merge local ids of facA and facB without using division information,
-  // i.e. append local ids of `facB` to `facA` and insert local ids of `facA`
-  // to `facB` at start of its local ids.
-  unsigned initLocals = facA.getNumLocalIds();
-  insertLocalId(facA.getNumLocalIds(), facB.getNumLocalIds());
-  facB.insertLocalId(0, initLocals);
-
-  // Get division representations from each FAC.
-  std::vector<SmallVector<int64_t, 8>> divsA, divsB;
-  SmallVector<unsigned, 4> denomsA, denomsB;
-  facA.getLocalReprs(divsA, denomsA);
-  facB.getLocalReprs(divsB, denomsB);
-
-  // Copy division information for facB into `divsA` and `denomsA`, so that
-  // these have the combined division information of both FACs. Since newly
-  // added local variables in facA and facB have no constraints, they will not
-  // have any division representation.
-  std::copy(divsB.begin() + initLocals, divsB.end(),
-            divsA.begin() + initLocals);
-  std::copy(denomsB.begin() + initLocals, denomsB.end(),
-            denomsA.begin() + initLocals);
-
-  // Find and merge duplicate divisions.
-  // TODO: Add division normalization to support divisions that differ by
-  // a constant.
-  // TODO: Add division ordering such that a division representation for local
-  // identifier at position `i` only depends on local identifiers at position <
-  // `i`. This would make sure that all divisions depending on other local
-  // variables that can be merged, are merged.
-  unsigned localOffset = getIdKindOffset(IdKind::Local);
-  for (unsigned i = 0; i < divsA.size(); ++i) {
-    // Check if a division representation exists for the `i^th` local id.
-    if (denomsA[i] == 0)
-      continue;
-    // Check if a division exists which is a duplicate of the division at `i`.
-    for (unsigned j = i + 1; j < divsA.size(); ++j) {
-      // Check if a division representation exists for the `j^th` local id.
-      if (denomsA[j] == 0)
-        continue;
-      // Check if the denominators match.
-      if (denomsA[i] != denomsA[j])
-        continue;
-      // Check if the representations are equal.
-      if (divsA[i] != divsA[j])
-        continue;
-
-      // Merge divisions at position `j` into division at position `i`.
-      eliminateRedundantLocalId(facA, i, j);
-      eliminateRedundantLocalId(facB, i, j);
-      for (unsigned k = 0, g = divsA.size(); k < g; ++k) {
-        SmallVector<int64_t, 8> &div = divsA[k];
-        if (denomsA[k] != 0) {
-          div[localOffset + i] += div[localOffset + j];
-          div.erase(div.begin() + localOffset + j);
-        }
-      }
-
-      divsA.erase(divsA.begin() + j);
-      denomsA.erase(denomsA.begin() + j);
-      // Since `j` can never be zero, we do not need to worry about overflows.
-      --j;
-    }
-  }
-}
-
-/// Removes local variables using equalities. Each equality is checked if it
-/// can be reduced to the form: `e = affine-expr`, where `e` is a local
-/// variable and `affine-expr` is an affine expression not containing `e`.
-/// If an equality satisfies this form, the local variable is replaced in
-/// each constraint and then removed. The equality used to replace this local
-/// variable is also removed.
-void FlatAffineConstraints::removeRedundantLocalVars() {
-  // Normalize the equality constraints to reduce coefficients of local
-  // variables to 1 wherever possible.
-  for (unsigned i = 0, e = getNumEqualities(); i < e; ++i)
-    normalizeConstraintByGCD</*isEq=*/true>(this, i);
-
-  while (true) {
-    unsigned i, e, j, f;
-    for (i = 0, e = getNumEqualities(); i < e; ++i) {
-      // Find a local variable to eliminate using ith equality.
-      for (j = getNumDimAndSymbolIds(), f = getNumIds(); j < f; ++j)
-        if (std::abs(atEq(i, j)) == 1)
-          break;
-
-      // Local variable can be eliminated using ith equality.
-      if (j < f)
-        break;
-    }
-
-    // No equality can be used to eliminate a local variable.
-    if (i == e)
-      break;
-
-    // Use the ith equality to simplify other equalities. If any changes
-    // are made to an equality constraint, it is normalized by GCD.
-    for (unsigned k = 0, t = getNumEqualities(); k < t; ++k) {
-      if (atEq(k, j) != 0) {
-        eliminateFromConstraint(this, k, i, j, j, /*isEq=*/true);
-        normalizeConstraintByGCD</*isEq=*/true>(this, k);
-      }
-    }
-
-    // Use the ith equality to simplify inequalities.
-    for (unsigned k = 0, t = getNumInequalities(); k < t; ++k)
-      eliminateFromConstraint(this, k, i, j, j, /*isEq=*/false);
-
-    // Remove the ith equality and the found local variable.
-    removeId(j);
-    removeEquality(i);
-  }
-}
-
-void FlatAffineConstraints::convertDimToLocal(unsigned dimStart,
-                                              unsigned dimLimit) {
-  assert(dimLimit <= getNumDimIds() && "Invalid dim pos range");
-
-  if (dimStart >= dimLimit)
-    return;
-
-  // Append new local variables corresponding to the dimensions to be converted.
-  unsigned convertCount = dimLimit - dimStart;
-  unsigned newLocalIdStart = getNumIds();
-  appendLocalId(convertCount);
-
-  // Swap the new local variables with dimensions.
-  for (unsigned i = 0; i < convertCount; ++i)
-    swapId(i + dimStart, i + newLocalIdStart);
-
-  // Remove dimensions converted to local variables.
-  removeIdRange(dimStart, dimLimit);
-}
-
 std::pair<AffineMap, AffineMap> FlatAffineConstraints::getLowerAndUpperBound(
     unsigned pos, unsigned offset, unsigned num, unsigned symStartPos,
     ArrayRef<AffineExpr> localExprs, MLIRContext *context) const {
   assert(pos + offset < getNumDimIds() && "invalid dim start pos");
   assert(symStartPos >= (pos + offset) && "invalid sym start pos");
   assert(getNumLocalIds() == localExprs.size() &&
          "incorrect local exprs count");
 
if (lbMap && lbMap.getNumResults() != 0 && ubMap &&
       auto loop = getForInductionVarOwner(values[i]);
       if (failed(this->addAffineForOpDomain(loop)))
         return failure();
     }
   }
   return success();
 }
 
-void FlatAffineConstraints::addBound(BoundType type, unsigned pos,
-                                     int64_t value) {
-  assert(pos < getNumCols());
-  if (type == BoundType::EQ) {
-    unsigned row = equalities.appendExtraRow();
-    equalities(row, pos) = 1;
-    equalities(row, getNumCols() - 1) = -value;
-  } else {
-    unsigned row = inequalities.appendExtraRow();
-    inequalities(row, pos) = type == BoundType::LB ? 1 : -1;
-    inequalities(row, getNumCols() - 1) =
-        type == BoundType::LB ? -value : value;
-  }
-}
-
-void FlatAffineConstraints::addBound(BoundType type, ArrayRef<int64_t> expr,
-                                     int64_t value) {
-  assert(type != BoundType::EQ && "EQ not implemented");
-  assert(expr.size() == getNumCols());
-  unsigned row = inequalities.appendExtraRow();
-  for (unsigned i = 0, e = expr.size(); i < e; ++i)
-    inequalities(row, i) = type == BoundType::LB ? expr[i] : -expr[i];
-  inequalities(inequalities.getNumRows() - 1, getNumCols() - 1) +=
-      type == BoundType::LB ? -value : value;
-}
-
-/// Adds a new local identifier as the floordiv of an affine function of other
-/// identifiers, the coefficients of which are provided in 'dividend' and with
-/// respect to a positive constant 'divisor'. Two constraints are added to the
-/// system to capture equivalence with the floordiv.
-///      q = expr floordiv c    <=>   c*q <= expr <= c*q + c - 1.
-void FlatAffineConstraints::addLocalFloorDiv(ArrayRef<int64_t> dividend,
-                                             int64_t divisor) {
-  assert(dividend.size() == getNumCols() && "incorrect dividend size");
-  assert(divisor > 0 && "positive divisor expected");
-
-  appendLocalId();
-
-  // Add two constraints for this new identifier 'q'.
-  SmallVector<int64_t, 8> bound(dividend.size() + 1);
-
-  // dividend - q * divisor >= 0
-  std::copy(dividend.begin(), dividend.begin() + dividend.size() - 1,
-            bound.begin());
-  bound.back() = dividend.back();
-  bound[getNumIds() - 1] = -divisor;
-  addInequality(bound);
-
-  // -dividend +qdivisor * q + divisor - 1 >= 0
-  std::transform(bound.begin(), bound.end(), bound.begin(),
-                 std::negate<int64_t>());
-  bound[bound.size() - 1] += divisor - 1;
-  addInequality(bound);
-}
-
 bool FlatAffineValueConstraints::findId(Value val, unsigned *pos) const {
   unsigned i = 0;
   for (const auto &mayBeId : values) {
     if (mayBeId.hasValue() && mayBeId.getValue() == val) {
       *pos = i;
       return true;
     }
     i++;
   }
   return false;
 }
 
 bool FlatAffineValueConstraints::containsId(Value val) const {
   return llvm::any_of(values, [&](const Optional<Value> &mayBeId) {
     return mayBeId.hasValue() && mayBeId.getValue() == val;
   });
 }
 
 void FlatAffineValueConstraints::swapId(unsigned posA, unsigned posB) {
   FlatAffineConstraints::swapId(posA, posB);
   std::swap(values[posA], values[posB]);
 }
 
-void FlatAffineConstraints::setDimSymbolSeparation(unsigned newSymbolCount) {
-  assert(newSymbolCount <= numDims + numSymbols &&
-         "invalid separation position");
-  numDims = numDims + numSymbols - newSymbolCount;
-  numSymbols = newSymbolCount;
-}
-
 void FlatAffineValueConstraints::addBound(BoundType type, Value val,
                                           int64_t value) {
   unsigned pos;
   if (!findId(val, &pos))
     // This is a pre-condition for this method.
     assert(0 && "id not found");
   addBound(type, pos, value);
 }
 
-/// Finds an equality that equates the specified identifier to a constant.
-/// Returns the position of the equality row. If 'symbolic' is set to true,
-/// symbols are also treated like a constant, i.e., an affine function of the
-/// symbols is also treated like a constant. Returns -1 if such an equality
-/// could not be found.
-static int findEqualityToConstant(const FlatAffineConstraints &cst,
-                                  unsigned pos, bool symbolic = false) {
-  assert(pos < cst.getNumIds() && "invalid position");
-  for (unsigned r = 0, e = cst.getNumEqualities(); r < e; r++) {
-    int64_t v = cst.atEq(r, pos);
-    if (v * v != 1)
-      continue;
-    unsigned c;
-    unsigned f = symbolic ? cst.getNumDimIds() : cst.getNumIds();
-    // This checks for zeros in all positions other than 'pos' in [0, f)
-    for (c = 0; c < f; c++) {
-      if (c == pos)
-        continue;
-      if (cst.atEq(r, c) != 0) {
-        // Dependent on another identifier.
-        break;
-      }
-    }
-    if (c == f)
-      // Equality is free of other identifiers.
-      return r;
-  }
-  return -1;
-}
-
-void FlatAffineConstraints::setAndEliminate(unsigned pos,
-                                            ArrayRef<int64_t> values) {
-  if (values.empty())
-    return;
-  assert(pos + values.size() <= getNumIds() &&
-         "invalid position or too many values");
-  // Setting x_j = p in sum_i a_i x_i + c is equivalent to adding p*a_j to the
-  // constant term and removing the id x_j. We do this for all the ids
-  // pos, pos + 1, ... pos + values.size() - 1.
-  for (unsigned r = 0, e = getNumInequalities(); r < e; r++)
-    for (unsigned i = 0, numVals = values.size(); i < numVals; ++i)
-      atIneq(r, getNumCols() - 1) += atIneq(r, pos + i) * values[i];
-  for (unsigned r = 0, e = getNumEqualities(); r < e; r++)
-    for (unsigned i = 0, numVals = values.size(); i < numVals; ++i)
-      atEq(r, getNumCols() - 1) += atEq(r, pos + i) * values[i];
-  removeIdRange(pos, pos + values.size());
-}
-
-LogicalResult FlatAffineConstraints::constantFoldId(unsigned pos) {
-  assert(pos < getNumIds() && "invalid position");
-  int rowIdx;
-  if ((rowIdx = findEqualityToConstant(*this, pos)) == -1)
-    return failure();
-
-  // atEq(rowIdx, pos) is either -1 or 1.
-  assert(atEq(rowIdx, pos) * atEq(rowIdx, pos) == 1);
-  int64_t constVal = -atEq(rowIdx, getNumCols() - 1) / atEq(rowIdx, pos);
-  setAndEliminate(pos, constVal);
-  return success();
-}
-
-void FlatAffineConstraints::constantFoldIdRange(unsigned pos, unsigned num) {
-  for (unsigned s = pos, t = pos, e = pos + num; s < e; s++) {
-    if (failed(constantFoldId(t)))
-      t++;
-  }
-}
-
-/// Returns a non-negative constant bound on the extent (upper bound - lower
-/// bound) of the specified identifier if it is found to be a constant; returns
-/// None if it's not a constant. This methods treats symbolic identifiers
-/// specially, i.e., it looks for constant differences between affine
-/// expressions involving only the symbolic identifiers. See comments at
-/// function definition for example. 'lb', if provided, is set to the lower
-/// bound associated with the constant difference. Note that 'lb' is purely
-/// symbolic and thus will contain the coefficients of the symbolic identifiers
-/// and the constant coefficient.
-//  Egs: 0 <= i <= 15, return 16.
-//       s0 + 2 <= i <= s0 + 17, returns 16. (s0 has to be a symbol)
-//       s0 + s1 + 16 <= d0 <= s0 + s1 + 31, returns 16.
-//       s0 - 7 <= 8*j <= s0 returns 1 with lb = s0, lbDivisor = 8 (since lb =
-//       ceil(s0 - 7 / 8) = floor(s0 / 8)).
-Optional<int64_t> FlatAffineConstraints::getConstantBoundOnDimSize(
-    unsigned pos, SmallVectorImpl<int64_t> *lb, int64_t *boundFloorDivisor,
-    SmallVectorImpl<int64_t> *ub, unsigned *minLbPos,
-    unsigned *minUbPos) const {
-  assert(pos < getNumDimIds() && "Invalid identifier position");
-
-  // Find an equality for 'pos'^th identifier that equates it to some function
-  // of the symbolic identifiers (+ constant).
-  int eqPos = findEqualityToConstant(*this, pos, /*symbolic=*/true);
-  if (eqPos != -1) {
-    auto eq = getEquality(eqPos);
-    // If the equality involves a local var, punt for now.
-    // TODO: this can be handled in the future by using the explicit
-    // representation of the local vars.
-    if (!std::all_of(eq.begin() + getNumDimAndSymbolIds(), eq.end() - 1,
-                     [](int64_t coeff) { return coeff == 0; }))
-      return None;
-
-    // This identifier can only take a single value.
-    if (lb) {
-      // Set lb to that symbolic value.
-      lb->resize(getNumSymbolIds() + 1);
-      if (ub)
-        ub->resize(getNumSymbolIds() + 1);
-      for (unsigned c = 0, f = getNumSymbolIds() + 1; c < f; c++) {
-        int64_t v = atEq(eqPos, pos);
-        // atEq(eqRow, pos) is either -1 or 1.
-        assert(v * v == 1);
-        (*lb)[c] = v < 0 ? atEq(eqPos, getNumDimIds() + c) / -v
-                         : -atEq(eqPos, getNumDimIds() + c) / v;
-        // Since this is an equality, ub = lb.
-        if (ub)
-          (*ub)[c] = (*lb)[c];
-      }
-      assert(boundFloorDivisor &&
-             "both lb and divisor or none should be provided");
-      *boundFloorDivisor = 1;
-    }
-    if (minLbPos)
-      *minLbPos = eqPos;
-    if (minUbPos)
-      *minUbPos = eqPos;
-    return 1;
-  }
-
-  // Check if the identifier appears at all in any of the inequalities.
-  unsigned r, e;
-  for (r = 0, e = getNumInequalities(); r < e; r++) {
-    if (atIneq(r, pos) != 0)
-      break;
-  }
-  if (r == e)
-    // If it doesn't, there isn't a bound on it.
-    return None;
-
-  // Positions of constraints that are lower/upper bounds on the variable.
-  SmallVector<unsigned, 4> lbIndices, ubIndices;
-
-  // Gather all symbolic lower bounds and upper bounds of the variable, i.e.,
-  // the bounds can only involve symbolic (and local) identifiers. Since the
-  // canonical form c_1*x_1 + c_2*x_2 + ... + c_0 >= 0, a constraint is a lower
-  // bound for x_i if c_i >= 1, and an upper bound if c_i <= -1.
-  getLowerAndUpperBoundIndices(pos, &lbIndices, &ubIndices,
-                               /*eqIndices=*/nullptr, /*offset=*/0,
-                               /*num=*/getNumDimIds());
-
-  Optional<int64_t> minDiff = None;
-  unsigned minLbPosition = 0, minUbPosition = 0;
-  for (auto ubPos : ubIndices) {
-    for (auto lbPos : lbIndices) {
-      // Look for a lower bound and an upper bound that only differ by a
-      // constant, i.e., pairs of the form  0 <= c_pos - f(c_i's) <= diffConst.
-      // For example, if ii is the pos^th variable, we are looking for
-      // constraints like ii >= i, ii <= ii + 50, 50 being the difference. The
-      // minimum among all such constant differences is kept since that's the
-      // constant bounding the extent of the pos^th variable.
-      unsigned j, e;
-      for (j = 0, e = getNumCols() - 1; j < e; j++)
-        if (atIneq(ubPos, j) != -atIneq(lbPos, j)) {
-          break;
-        }
-      if (j < getNumCols() - 1)
-        continue;
-      int64_t diff = ceilDiv(atIneq(ubPos, getNumCols() - 1) +
-                                 atIneq(lbPos, getNumCols() - 1) + 1,
-                             atIneq(lbPos, pos));
-      // This bound is non-negative by definition.
-      diff = std::max<int64_t>(diff, 0);
-      if (minDiff == None || diff < minDiff) {
-        minDiff = diff;
-        minLbPosition = lbPos;
-        minUbPosition = ubPos;
-      }
-    }
-  }
-  if (lb && minDiff.hasValue()) {
-    // Set lb to the symbolic lower bound.
-    lb->resize(getNumSymbolIds() + 1);
-    if (ub)
-      ub->resize(getNumSymbolIds() + 1);
-    // The lower bound is the ceildiv of the lb constraint over the coefficient
-    // of the variable at 'pos'. We express the ceildiv equivalently as a floor
-    // for uniformity. For eg., if the lower bound constraint was: 32*d0 - N +
-    // 31 >= 0, the lower bound for d0 is ceil(N - 31, 32), i.e., floor(N, 32).
-    *boundFloorDivisor = atIneq(minLbPosition, pos);
-    assert(*boundFloorDivisor == -atIneq(minUbPosition, pos));
-    for (unsigned c = 0, e = getNumSymbolIds() + 1; c < e; c++) {
-      (*lb)[c] = -atIneq(minLbPosition, getNumDimIds() + c);
-    }
-    if (ub) {
-      for (unsigned c = 0, e = getNumSymbolIds() + 1; c < e; c++)
-        (*ub)[c] = atIneq(minUbPosition, getNumDimIds() + c);
-    }
-    // The lower bound leads to a ceildiv while the upper bound is a floordiv
-    // whenever the coefficient at pos != 1. ceildiv (val / d) = floordiv (val +
-    // d - 1 / d); hence, the addition of 'atIneq(minLbPosition, pos) - 1' to
-    // the constant term for the lower bound.
-    (*lb)[getNumSymbolIds()] += atIneq(minLbPosition, pos) - 1;
-  }
-  if (minLbPos)
-    *minLbPos = minLbPosition;
-  if (minUbPos)
-    *minUbPos = minUbPosition;
-  return minDiff;
-}
-
-template <bool isLower>
-Optional<int64_t>
-FlatAffineConstraints::computeConstantLowerOrUpperBound(unsigned pos) {
-  assert(pos < getNumIds() && "invalid position");
-  // Project to 'pos'.
-  projectOut(0, pos);
-  projectOut(1, getNumIds() - 1);
-  // Check if there's an equality equating the '0'^th identifier to a constant.
-  int eqRowIdx = findEqualityToConstant(*this, 0, /*symbolic=*/false);
-  if (eqRowIdx != -1)
-    // atEq(rowIdx, 0) is either -1 or 1.
-    return -atEq(eqRowIdx, getNumCols() - 1) / atEq(eqRowIdx, 0);
-
-  // Check if the identifier appears at all in any of the inequalities.
-  unsigned r, e;
-  for (r = 0, e = getNumInequalities(); r < e; r++) {
-    if (atIneq(r, 0) != 0)
-      break;
-  }
-  if (r == e)
-    // If it doesn't, there isn't a bound on it.
-    return None;
-
-  Optional<int64_t> minOrMaxConst = None;
-
-  // Take the max across all const lower bounds (or min across all constant
-  // upper bounds).
-  for (unsigned r = 0, e = getNumInequalities(); r < e; r++) {
-    if (isLower) {
-      if (atIneq(r, 0) <= 0)
-        // Not a lower bound.
-        continue;
-    } else if (atIneq(r, 0) >= 0) {
-      // Not an upper bound.
-      continue;
-    }
-    unsigned c, f;
-    for (c = 0, f = getNumCols() - 1; c < f; c++)
-      if (c != 0 && atIneq(r, c) != 0)
-        break;
-    if (c < getNumCols() - 1)
-      // Not a constant bound.
-      continue;
-
-    int64_t boundConst =
-        isLower ? mlir::ceilDiv(-atIneq(r, getNumCols() - 1), atIneq(r, 0))
-                : mlir::floorDiv(atIneq(r, getNumCols() - 1), -atIneq(r, 0));
-    if (isLower) {
-      if (minOrMaxConst == None || boundConst > minOrMaxConst)
-        minOrMaxConst = boundConst;
-    } else {
-      if (minOrMaxConst == None || boundConst < minOrMaxConst)
-        minOrMaxConst = boundConst;
-    }
-  }
-  return minOrMaxConst;
-}
-
-Optional<int64_t> FlatAffineConstraints::getConstantBound(BoundType type,
-                                                          unsigned pos) const {
-  assert(type != BoundType::EQ && "EQ not implemented");
-  FlatAffineConstraints tmpCst(*this);
-  if (type == BoundType::LB)
-    return tmpCst.computeConstantLowerOrUpperBound</*isLower=*/true>(pos);
-  return tmpCst.computeConstantLowerOrUpperBound</*isLower=*/false>(pos);
-}
-
-// A simple (naive and conservative) check for hyper-rectangularity.
-bool FlatAffineConstraints::isHyperRectangular(unsigned pos,
-                                               unsigned num) const {
-  assert(pos < getNumCols() - 1);
-  // Check for two non-zero coefficients in the range [pos, pos + sum).
-  for (unsigned r = 0, e = getNumInequalities(); r < e; r++) {
-    unsigned sum = 0;
-    for (unsigned c = pos; c < pos + num; c++) {
-      if (atIneq(r, c) != 0)
-        sum++;
-    }
-    if (sum > 1)
-      return false;
-  }
-  for (unsigned r = 0, e = getNumEqualities(); r < e; r++) {
-    unsigned sum = 0;
-    for (unsigned c = pos; c < pos + num; c++) {
-      if (atEq(r, c) != 0)
-        sum++;
-    }
-    if (sum > 1)
-      return false;
-  }
-  return true;
-}
-
-void FlatAffineConstraints::print(raw_ostream &os) const {
-  assert(hasConsistentState());
-  os << "\nConstraints (" << getNumDimIds() << " dims, " << getNumSymbolIds()
-     << " symbols, " << getNumLocalIds() << " locals), (" << getNumConstraints()
-     << " constraints)\n";
+void FlatAffineConstraints::printSpace(raw_ostream &os) const {
+  IntegerPolyhedron::print(os);
   os << "(";
   for (unsigned i = 0, e = getNumIds(); i < e; i++) {
     if (auto *valueCstr = dyn_cast<const FlatAffineValueConstraints>(this)) {
       if (valueCstr->hasValue(i))
         os << "Value ";
       else
         os << "None ";
     } else {
       os << "None ";
     }
   }
   os << " const)\n";
-  for (unsigned i = 0, e = getNumEqualities(); i < e; ++i) {
-    for (unsigned j = 0, f = getNumCols(); j < f; ++j) {
-      os << atEq(i, j) << " ";
-    }
-    os << "= 0\n";
-  }
-  for (unsigned i = 0, e = getNumInequalities(); i < e; ++i) {
-    for (unsigned j = 0, f = getNumCols(); j < f; ++j) {
-      os << atIneq(i, j) << " ";
-    }
-    os << ">= 0\n";
-  }
-  os << '\n';
-}
-
-void FlatAffineConstraints::dump() const { print(llvm::errs()); }
-
-/// Removes duplicate constraints, trivially true constraints, and constraints
-/// that can be detected as redundant as a result of differing only in their
-/// constant term part. A constraint of the form <non-negative constant> >= 0 is
-/// considered trivially true.
-//  Uses a DenseSet to hash and detect duplicates followed by a linear scan to
-//  remove duplicates in place.
-void FlatAffineConstraints::removeTrivialRedundancy() {
-  gcdTightenInequalities();
-  normalizeConstraintsByGCD();
-
-  // A map used to detect redundancy stemming from constraints that only differ
-  // in their constant term. The value stored is <row position, const term>
-  // for a given row.
-  SmallDenseMap<ArrayRef<int64_t>, std::pair<unsigned, int64_t>>
-      rowsWithoutConstTerm;
-  // To unique rows.
-  SmallDenseSet<ArrayRef<int64_t>, 8> rowSet;
-
-  // Check if constraint is of the form <non-negative-constant> >= 0.
-  auto isTriviallyValid = [&](unsigned r) -> bool {
-    for (unsigned c = 0, e = getNumCols() - 1; c < e; c++) {
-      if (atIneq(r, c) != 0)
-        return false;
-    }
-    return atIneq(r, getNumCols() - 1) >= 0;
-  };
-
-  // Detect and mark redundant constraints.
-  SmallVector<bool, 256> redunIneq(getNumInequalities(), false);
-  for (unsigned r = 0, e = getNumInequalities(); r < e; r++) {
-    int64_t *rowStart = &inequalities(r, 0);
-    auto row = ArrayRef<int64_t>(rowStart, getNumCols());
-    if (isTriviallyValid(r) || !rowSet.insert(row).second) {
-      redunIneq[r] = true;
-      continue;
-    }
-
-    // Among constraints that only differ in the constant term part, mark
-    // everything other than the one with the smallest constant term redundant.
-    // (eg: among i - 16j - 5 >= 0, i - 16j - 1 >=0, i - 16j - 7 >= 0, the
-    // former two are redundant).
-    int64_t constTerm = atIneq(r, getNumCols() - 1);
-    auto rowWithoutConstTerm = ArrayRef<int64_t>(rowStart, getNumCols() - 1);
-    const auto &ret =
-        rowsWithoutConstTerm.insert({rowWithoutConstTerm, {r, constTerm}});
-    if (!ret.second) {
-      // Check if the other constraint has a higher constant term.
-      auto &val = ret.first->second;
-      if (val.second > constTerm) {
-        // The stored row is redundant. Mark it so, and update with this one.
-        redunIneq[val.first] = true;
-        val = {r, constTerm};
-      } else {
-        // The one stored makes this one redundant.
-        redunIneq[r] = true;
-      }
-    }
-  }
-
-  // Scan to get rid of all rows marked redundant, in-place.
-  unsigned pos = 0;
-  for (unsigned r = 0, e = getNumInequalities(); r < e; r++)
-    if (!redunIneq[r])
-      inequalities.copyRow(r, pos++);
-
-  inequalities.resizeVertically(pos);
 
-  // TODO: consider doing this for equalities as well, but probably not worth
-  // the savings.
-}
-
-void FlatAffineConstraints::clearAndCopyFrom(
-    const FlatAffineConstraints &other) {
-  if (auto *otherValueSet = dyn_cast<const FlatAffineValueConstraints>(&other))
+void FlatAffineConstraints::clearAndCopyFrom(const IntegerPolyhedron &other) {
+  const FlatAffineConstraints *otherFAC =
+      static_cast<const FlatAffineConstraints *>(&other);
+  if (auto *otherValueSet =
+          dyn_cast<const FlatAffineValueConstraints>(otherFAC))
     assert(!otherValueSet->hasValues() &&
            "cannot copy associated Values into FlatAffineConstraints");
   // Note: Assigment operator does not vtable pointer, so kind does not change.
-  *this = other;
+  *static_cast<FlatAffineConstraints *>(this) = *otherFAC;
 }
 
 void FlatAffineValueConstraints::clearAndCopyFrom(
-    const FlatAffineConstraints &other) {
+    const IntegerPolyhedron &other) {
+  const FlatAffineConstraints *otherFAC =
+      static_cast<const FlatAffineConstraints *>(&other);
   if (auto *otherValueSet =
-          dyn_cast<const FlatAffineValueConstraints>(&other)) {
+          dyn_cast<const FlatAffineValueConstraints>(otherFAC)) {
     *this = *otherValueSet;
   } else {
-    *static_cast<FlatAffineConstraints *>(this) = other;
+    *static_cast<FlatAffineConstraints *>(this) = *otherFAC;
     values.clear();
     values.resize(numIds, None);
   }
 }
 
-static std::pair<unsigned, unsigned>
-getNewNumDimsSymbols(unsigned pos, const FlatAffineConstraints &cst) {
-  unsigned numDims = cst.getNumDimIds();
-  unsigned numSymbols = cst.getNumSymbolIds();
-  unsigned newNumDims, newNumSymbols;
-  if (pos < numDims) {
-    newNumDims = numDims - 1;
-    newNumSymbols = numSymbols;
-  } else if (pos < numDims + numSymbols) {
-    assert(numSymbols >= 1);
-    newNumDims = numDims;
-    newNumSymbols = numSymbols - 1;
-  } else {
-    newNumDims = numDims;
-    newNumSymbols = numSymbols;
-  }
-  return {newNumDims, newNumSymbols};
-}
-
-#undef DEBUG_TYPE
-#define DEBUG_TYPE "fm"
-
-/// Eliminates identifier at the specified position using Fourier-Motzkin
-/// variable elimination. This technique is exact for rational spaces but
-/// conservative (in "rare" cases) for integer spaces. The operation corresponds
-/// to a projection operation yielding the (convex) set of integer points
-/// contained in the rational shadow of the set. An emptiness test that relies
-/// on this method will guarantee emptiness, i.e., it disproves the existence of
-/// a solution if it says it's empty.
-/// If a non-null isResultIntegerExact is passed, it is set to true if the
-/// result is also integer exact. If it's set to false, the obtained solution
-/// *may* not be exact, i.e., it may contain integer points that do not have an
-/// integer pre-image in the original set.
-///
-/// Eg:
-/// j >= 0, j <= i + 1
-/// i >= 0, i <= N + 1
-/// Eliminating i yields,
-///   j >= 0, 0 <= N + 1, j - 1 <= N + 1
-///
-/// If darkShadow = true, this method computes the dark shadow on elimination;
-/// the dark shadow is a convex integer subset of the exact integer shadow. A
-/// non-empty dark shadow proves the existence of an integer solution. The
-/// elimination in such a case could however be an under-approximation, and thus
-/// should not be used for scanning sets or used by itself for dependence
-/// checking.
-///
-/// Eg: 2-d set, * represents grid points, 'o' represents a point in the set.
-///            ^
-///            |
-///            | * * * * o o
-///         i  | * * o o o o
-///            | o * * * * *
-///            --------------->
-///                 j ->
-///
-/// Eliminating i from this system (projecting on the j dimension):
-/// rational shadow / integer light shadow:  1 <= j <= 6
-/// dark shadow:                             3 <= j <= 6
-/// exact integer shadow:                    j = 1 \union  3 <= j <= 6
-/// holes/splinters:                         j = 2
-///
-/// darkShadow = false, isResultIntegerExact = nullptr are default values.
-// TODO: a slight modification to yield dark shadow version of FM (tightened),
-// which can prove the existence of a solution if there is one.
-void FlatAffineConstraints::fourierMotzkinEliminate(
-    unsigned pos, bool darkShadow, bool *isResultIntegerExact) {
-  LLVM_DEBUG(llvm::dbgs() << "FM input (eliminate pos " << pos << "):\n");
-  LLVM_DEBUG(dump());
-  assert(pos < getNumIds() && "invalid position");
-  assert(hasConsistentState());
-
-  // Check if this identifier can be eliminated through a substitution.
-  for (unsigned r = 0, e = getNumEqualities(); r < e; r++) {
-    if (atEq(r, pos) != 0) {
-      // Use Gaussian elimination here (since we have an equality).
-      LogicalResult ret = gaussianEliminateId(pos);
-      (void)ret;
-      assert(succeeded(ret) && "Gaussian elimination guaranteed to succeed");
-      LLVM_DEBUG(llvm::dbgs() << "FM output (through Gaussian elimination):\n");
-      LLVM_DEBUG(dump());
-      return;
-    }
-  }
-
-  // A fast linear time tightening.
-  gcdTightenInequalities();
-
-  // Check if the identifier appears at all in any of the inequalities.
-  unsigned r, e;
-  for (r = 0, e = getNumInequalities(); r < e; r++) {
-    if (atIneq(r, pos) != 0)
-      break;
-  }
-  if (r == getNumInequalities()) {
-    // If it doesn't appear, just remove the column and return.
-    // TODO: refactor removeColumns to use it from here.
-    removeId(pos);
-    LLVM_DEBUG(llvm::dbgs() << "FM output:\n");
-    LLVM_DEBUG(dump());
-    return;
-  }
-
-  // Positions of constraints that are lower bounds on the variable.
-  SmallVector<unsigned, 4> lbIndices;
-  // Positions of constraints that are lower bounds on the variable.
-  SmallVector<unsigned, 4> ubIndices;
-  // Positions of constraints that do not involve the variable.
-  std::vector<unsigned> nbIndices;
-  nbIndices.reserve(getNumInequalities());
-
-  // Gather all lower bounds and upper bounds of the variable. Since the
-  // canonical form c_1*x_1 + c_2*x_2 + ... + c_0 >= 0, a constraint is a lower
-  // bound for x_i if c_i >= 1, and an upper bound if c_i <= -1.
-  for (unsigned r = 0, e = getNumInequalities(); r < e; r++) {
-    if (atIneq(r, pos) == 0) {
-      // Id does not appear in bound.
-      nbIndices.push_back(r);
-    } else if (atIneq(r, pos) >= 1) {
-      // Lower bound.
-      lbIndices.push_back(r);
-    } else {
-      // Upper bound.
-      ubIndices.push_back(r);
-    }
-  }
-
-  // Set the number of dimensions, symbols in the resulting system.
-  const auto &dimsSymbols = getNewNumDimsSymbols(pos, *this);
-  unsigned newNumDims = dimsSymbols.first;
-  unsigned newNumSymbols = dimsSymbols.second;
-
-  /// Create the new system which has one identifier less.
-  FlatAffineConstraints newFac(
-      lbIndices.size() * ubIndices.size() + nbIndices.size(),
-      getNumEqualities(), getNumCols() - 1, newNumDims, newNumSymbols,
-      /*numLocals=*/getNumIds() - 1 - newNumDims - newNumSymbols);
-
-  // This will be used to check if the elimination was integer exact.
-  unsigned lcmProducts = 1;
-
-  // Let x be the variable we are eliminating.
-  // For each lower bound, lb <= c_l*x, and each upper bound c_u*x <= ub, (note
-  // that c_l, c_u >= 1) we have:
-  // lb*lcm(c_l, c_u)/c_l <= lcm(c_l, c_u)*x <= ub*lcm(c_l, c_u)/c_u
-  // We thus generate a constraint:
-  // lcm(c_l, c_u)/c_l*lb <= lcm(c_l, c_u)/c_u*ub.
-  // Note if c_l = c_u = 1, all integer points captured by the resulting
-  // constraint correspond to integer points in the original system (i.e., they
-  // have integer pre-images). Hence, if the lcm's are all 1, the elimination is
-  // integer exact.
-  for (auto ubPos : ubIndices) {
-    for (auto lbPos : lbIndices) {
-      SmallVector<int64_t, 4> ineq;
-      ineq.reserve(newFac.getNumCols());
-      int64_t lbCoeff = atIneq(lbPos, pos);
-      // Note that in the comments above, ubCoeff is the negation of the
-      // coefficient in the canonical form as the view taken here is that of the
-      // term being moved to the other size of '>='.
-      int64_t ubCoeff = -atIneq(ubPos, pos);
-      // TODO: refactor this loop to avoid all branches inside.
-      for (unsigned l = 0, e = getNumCols(); l < e; l++) {
-        if (l == pos)
-          continue;
-        assert(lbCoeff >= 1 && ubCoeff >= 1 && "bounds wrongly identified");
-        int64_t lcm = mlir::lcm(lbCoeff, ubCoeff);
-        ineq.push_back(atIneq(ubPos, l) * (lcm / ubCoeff) +
-                       atIneq(lbPos, l) * (lcm / lbCoeff));
-        lcmProducts *= lcm;
-      }
-      if (darkShadow) {
-        // The dark shadow is a convex subset of the exact integer shadow. If
-        // there is a point here, it proves the existence of a solution.
-        ineq[ineq.size() - 1] += lbCoeff * ubCoeff - lbCoeff - ubCoeff + 1;
-      }
-      // TODO: we need to have a way to add inequalities in-place in
-      // FlatAffineConstraints instead of creating and copying over.
-      newFac.addInequality(ineq);
-    }
-  }
-
-  LLVM_DEBUG(llvm::dbgs() << "FM isResultIntegerExact: " << (lcmProducts == 1)
-                          << "\n");
-  if (lcmProducts == 1 && isResultIntegerExact)
-    *isResultIntegerExact = true;
-
-  // Copy over the constraints not involving this variable.
-  for (auto nbPos : nbIndices) {
-    SmallVector<int64_t, 4> ineq;
-    ineq.reserve(getNumCols() - 1);
-    for (unsigned l = 0, e = getNumCols(); l < e; l++) {
-      if (l == pos)
-        continue;
-      ineq.push_back(atIneq(nbPos, l));
-    }
-    newFac.addInequality(ineq);
-  }
-
-  assert(newFac.getNumConstraints() ==
-         lbIndices.size() * ubIndices.size() + nbIndices.size());
-
-  // Copy over the equalities.
-  for (unsigned r = 0, e = getNumEqualities(); r < e; r++) {
-    SmallVector<int64_t, 4> eq;
-    eq.reserve(newFac.getNumCols());
-    for (unsigned l = 0, e = getNumCols(); l < e; l++) {
-      if (l == pos)
-        continue;
-      eq.push_back(atEq(r, l));
-    }
-    newFac.addEquality(eq);
-  }
-
-  // GCD tightening and normalization allows detection of more trivially
-  // redundant constraints.
-  newFac.gcdTightenInequalities();
-  newFac.normalizeConstraintsByGCD();
-  newFac.removeTrivialRedundancy();
-  clearAndCopyFrom(newFac);
-  LLVM_DEBUG(llvm::dbgs() << "FM output:\n");
-  LLVM_DEBUG(dump());
-}
-
-#undef DEBUG_TYPE
-#define DEBUG_TYPE "affine-structures"
-
 void FlatAffineValueConstraints::fourierMotzkinEliminate(
     unsigned pos, bool darkShadow, bool *isResultIntegerExact) {
   SmallVector<Optional<Value>, 8> newVals;
   newVals.reserve(numIds - 1);
   newVals.append(values.begin(), values.begin() + pos);
   newVals.append(values.begin() + pos + 1, values.end());
   // Note: Base implementation discards all associated Values.
   FlatAffineConstraints::fourierMotzkinEliminate(pos, darkShadow,
                                                  isResultIntegerExact);
   values = newVals;
   assert(values.size() == getNumIds());
 }
 
-void FlatAffineConstraints::projectOut(unsigned pos, unsigned num) {
-  if (num == 0)
-    return;
-
-  // 'pos' can be at most getNumCols() - 2 if num > 0.
-  assert((getNumCols() < 2 || pos <= getNumCols() - 2) && "invalid position");
-  assert(pos + num < getNumCols() && "invalid range");
-
-  // Eliminate as many identifiers as possible using Gaussian elimination.
-  unsigned currentPos = pos;
-  unsigned numToEliminate = num;
-  unsigned numGaussianEliminated = 0;
-
-  while (currentPos < getNumIds()) {
-    unsigned curNumEliminated =
-        gaussianEliminateIds(currentPos, currentPos + numToEliminate);
-    ++currentPos;
-    numToEliminate -= curNumEliminated + 1;
-    numGaussianEliminated += curNumEliminated;
-  }
-
-  // Eliminate the remaining using Fourier-Motzkin.
-  for (unsigned i = 0; i < num - numGaussianEliminated; i++) {
-    unsigned numToEliminate = num - numGaussianEliminated - i;
-    fourierMotzkinEliminate(
-        getBestIdToEliminate(*this, pos, pos + numToEliminate));
-  }
-
-  // Fast/trivial simplifications.
-  gcdTightenInequalities();
-  // Normalize constraints after tightening since the latter impacts this, but
-  // not the other way round.
-  normalizeConstraintsByGCD();
-}
-
 void FlatAffineValueConstraints::projectOut(Value val) {
   unsigned pos;
   bool ret = findId(val, &pos);
   assert(ret);
   (void)ret;
   fourierMotzkinEliminate(pos);
 }
 
-namespace {
-
-enum BoundCmpResult { Greater, Less, Equal, Unknown };
-
-/// Compares two affine bounds whose coefficients are provided in 'first' and
-/// 'second'. The last coefficient is the constant term.
-static BoundCmpResult compareBounds(ArrayRef<int64_t> a, ArrayRef<int64_t> b) {
-  assert(a.size() == b.size());
-
-  // For the bounds to be comparable, their corresponding identifier
-  // coefficients should be equal; the constant terms are then compared to
-  // determine less/greater/equal.
-
-  if (!std::equal(a.begin(), a.end() - 1, b.begin()))
-    return Unknown;
-
-  if (a.back() == b.back())
-    return Equal;
-
-  return a.back() < b.back() ? Less : Greater;
-}
-} // namespace
-
-// Returns constraints that are common to both A & B.
-static void getCommonConstraints(const FlatAffineConstraints &a,
-                                 const FlatAffineConstraints &b,
-                                 FlatAffineConstraints &c) {
-  c.reset(a.getNumDimIds(), a.getNumSymbolIds(), a.getNumLocalIds());
-  // a naive O(n^2) check should be enough here given the input sizes.
-  for (unsigned r = 0, e = a.getNumInequalities(); r < e; ++r) {
-    for (unsigned s = 0, f = b.getNumInequalities(); s < f; ++s) {
-      if (a.getInequality(r) == b.getInequality(s)) {
-        c.addInequality(a.getInequality(r));
-        break;
-      }
-    }
-  }
-  for (unsigned r = 0, e = a.getNumEqualities(); r < e; ++r) {
-    for (unsigned s = 0, f = b.getNumEqualities(); s < f; ++s) {
-      if (a.getEquality(r) == b.getEquality(s)) {
-        c.addEquality(a.getEquality(r));
-        break;
-      }
-    }
-  }
-}
-
-// Computes the bounding box with respect to 'other' by finding the min of the
-// lower bounds and the max of the upper bounds along each of the dimensions.
-LogicalResult
-FlatAffineConstraints::unionBoundingBox(const FlatAffineConstraints &otherCst) {
-  assert(otherCst.getNumDimIds() == numDims && "dims mismatch");
-  assert(otherCst.getNumLocalIds() == 0 && "local ids not supported here");
-  assert(getNumLocalIds() == 0 && "local ids not supported yet here");
-
-  // Get the constraints common to both systems; these will be added as is to
-  // the union.
-  FlatAffineConstraints commonCst;
-  getCommonConstraints(*this, otherCst, commonCst);
-
-  std::vector<SmallVector<int64_t, 8>> boundingLbs;
-  std::vector<SmallVector<int64_t, 8>> boundingUbs;
-  boundingLbs.reserve(2 * getNumDimIds());
-  boundingUbs.reserve(2 * getNumDimIds());
-
-  // To hold lower and upper bounds for each dimension.
-  SmallVector<int64_t, 4> lb, otherLb, ub, otherUb;
-  // To compute min of lower bounds and max of upper bounds for each dimension.
-  SmallVector<int64_t, 4> minLb(getNumSymbolIds() + 1);
-  SmallVector<int64_t, 4> maxUb(getNumSymbolIds() + 1);
-  // To compute final new lower and upper bounds for the union.
-  SmallVector<int64_t, 8> newLb(getNumCols()), newUb(getNumCols());
-
-  int64_t lbFloorDivisor, otherLbFloorDivisor;
-  for (unsigned d = 0, e = getNumDimIds(); d < e; ++d) {
-    auto extent = getConstantBoundOnDimSize(d, &lb, &lbFloorDivisor, &ub);
-    if (!extent.hasValue())
-      // TODO: symbolic extents when necessary.
-      // TODO: handle union if a dimension is unbounded.
-      return failure();
-
-    auto otherExtent = otherCst.getConstantBoundOnDimSize(
-        d, &otherLb, &otherLbFloorDivisor, &otherUb);
-    if (!otherExtent.hasValue() || lbFloorDivisor != otherLbFloorDivisor)
-      // TODO: symbolic extents when necessary.
-      return failure();
-
-    assert(lbFloorDivisor > 0 && "divisor always expected to be positive");
-
-    auto res = compareBounds(lb, otherLb);
-    // Identify min.
-    if (res == BoundCmpResult::Less || res == BoundCmpResult::Equal) {
-      minLb = lb;
-      // Since the divisor is for a floordiv, we need to convert to ceildiv,
-      // i.e., i >= expr floordiv div <=> i >= (expr - div + 1) ceildiv div <=>
-      // div * i >= expr - div + 1.
-      minLb.back() -= lbFloorDivisor - 1;
-    } else if (res == BoundCmpResult::Greater) {
-      minLb = otherLb;
-      minLb.back() -= otherLbFloorDivisor - 1;
-    } else {
-      // Uncomparable - check for constant lower/upper bounds.
-      auto constLb = getConstantBound(BoundType::LB, d);
-      auto constOtherLb = otherCst.getConstantBound(BoundType::LB, d);
-      if (!constLb.hasValue() || !constOtherLb.hasValue())
-        return failure();
-      std::fill(minLb.begin(), minLb.end(), 0);
-      minLb.back() = std::min(constLb.getValue(), constOtherLb.getValue());
-    }
-
-    // Do the same for ub's but max of upper bounds. Identify max.
-    auto uRes = compareBounds(ub, otherUb);
-    if (uRes == BoundCmpResult::Greater || uRes == BoundCmpResult::Equal) {
-      maxUb = ub;
-    } else if (uRes == BoundCmpResult::Less) {
-      maxUb = otherUb;
-    } else {
-      // Uncomparable - check for constant lower/upper bounds.
-      auto constUb = getConstantBound(BoundType::UB, d);
-      auto constOtherUb = otherCst.getConstantBound(BoundType::UB, d);
-      if (!constUb.hasValue() || !constOtherUb.hasValue())
-        return failure();
-      std::fill(maxUb.begin(), maxUb.end(), 0);
-      maxUb.back() = std::max(constUb.getValue(), constOtherUb.getValue());
-    }
-
-    std::fill(newLb.begin(), newLb.end(), 0);
-    std::fill(newUb.begin(), newUb.end(), 0);
-
-    // The divisor for lb, ub, otherLb, otherUb at this point is lbDivisor,
-    // and so it's the divisor for newLb and newUb as well.
-    newLb[d] = lbFloorDivisor;
-    newUb[d] = -lbFloorDivisor;
-    // Copy over the symbolic part + constant term.
-    std::copy(minLb.begin(), minLb.end(), newLb.begin() + getNumDimIds());
-    std::transform(newLb.begin() + getNumDimIds(), newLb.end(),
-                   newLb.begin() + getNumDimIds(), std::negate<int64_t>());
-    std::copy(maxUb.begin(), maxUb.end(), newUb.begin() + getNumDimIds());
-
-    boundingLbs.push_back(newLb);
-    boundingUbs.push_back(newUb);
-  }
-
-  // Clear all constraints and add the lower/upper bounds for the bounding box.
-  clearConstraints();
-  for (unsigned d = 0, e = getNumDimIds(); d < e; ++d) {
-    addInequality(boundingLbs[d]);
-    addInequality(boundingUbs[d]);
-  }
-
-  // Add the constraints that were common to both systems.
-  append(commonCst);
-  removeTrivialRedundancy();
-
-  // TODO: copy over pure symbolic constraints from this and 'other' over to the
-  // union (since the above are just the union along dimensions); we shouldn't
-  // be discarding any other constraints on the symbols.
-
-  return success();
-}
-
 LogicalResult FlatAffineValueConstraints::unionBoundingBox(
     const FlatAffineValueConstraints &otherCst) {
   assert(otherCst.getNumDimIds() == numDims && "dims mismatch");
   assert(otherCst.getMaybeValues()
              .slice(0, getNumDimIds())
              .equals(getMaybeValues().slice(0, getNumDimIds())) &&
          "dim values mismatch");
   assert(otherCst.getNumLocalIds() == 0 && "local ids not supported here");
void FlatAffineValueConstraints::getIneqAsAffineValueMap(
   SmallVector<Value, 4> operands;
   getValues(0, pos, &operands);
   SmallVector<Value, 4> trailingOperands;
   getValues(pos + 1, getNumDimAndSymbolIds(), &trailingOperands);
   operands.append(trailingOperands.begin(), trailingOperands.end());
   vmap.reset(AffineMap::get(numDims - 1, numSyms, boundExpr), operands);
 }
 
-bool FlatAffineConstraints::isColZero(unsigned pos) const {
-  unsigned rowPos;
-  return !findConstraintWithNonZeroAt(pos, /*isEq=*/false, &rowPos) &&
-         !findConstraintWithNonZeroAt(pos, /*isEq=*/true, &rowPos);
-}
-
 IntegerSet FlatAffineConstraints::getAsIntegerSet(MLIRContext *context) const {
   if (getNumConstraints() == 0)
     // Return universal set (always true): 0 == 0.
     return IntegerSet::get(getNumDimIds(), getNumSymbolIds(),
                            getAffineConstantExpr(/*constant=*/0, context),
                            /*eqFlags=*/true);
 
   // Construct local references.
for (unsigned i = 0, e = getNumEqualities(); i < e; ++i)
     exprs.push_back(getAffineExprFromFlatForm(getEquality(i), numDims, numSyms,
                                               localExprs, context));
   for (unsigned i = 0, e = getNumInequalities(); i < e; ++i)
     exprs.push_back(getAffineExprFromFlatForm(getInequality(i), numDims,
                                               numSyms, localExprs, context));
   return IntegerSet::get(numDims, numSyms, exprs, eqFlags);
 }
 
-/// Find positions of inequalities and equalities that do not have a coefficient
-/// for [pos, pos + num) identifiers.
-static void getIndependentConstraints(const FlatAffineConstraints &cst,
-                                      unsigned pos, unsigned num,
-                                      SmallVectorImpl<unsigned> &nbIneqIndices,
-                                      SmallVectorImpl<unsigned> &nbEqIndices) {
-  assert(pos < cst.getNumIds() && "invalid start position");
-  assert(pos + num <= cst.getNumIds() && "invalid limit");
-
-  for (unsigned r = 0, e = cst.getNumInequalities(); r < e; r++) {
-    // The bounds are to be independent of [offset, offset + num) columns.
-    unsigned c;
-    for (c = pos; c < pos + num; ++c) {
-      if (cst.atIneq(r, c) != 0)
-        break;
-    }
-    if (c == pos + num)
-      nbIneqIndices.push_back(r);
-  }
-
-  for (unsigned r = 0, e = cst.getNumEqualities(); r < e; r++) {
-    // The bounds are to be independent of [offset, offset + num) columns.
-    unsigned c;
-    for (c = pos; c < pos + num; ++c) {
-      if (cst.atEq(r, c) != 0)
-        break;
-    }
-    if (c == pos + num)
-      nbEqIndices.push_back(r);
-  }
-}
-
-void FlatAffineConstraints::removeIndependentConstraints(unsigned pos,
-                                                         unsigned num) {
-  assert(pos + num <= getNumIds() && "invalid range");
-
-  // Remove constraints that are independent of these identifiers.
-  SmallVector<unsigned, 4> nbIneqIndices, nbEqIndices;
-  getIndependentConstraints(*this, /*pos=*/0, num, nbIneqIndices, nbEqIndices);
-
-  // Iterate in reverse so that indices don't have to be updated.
-  // TODO: This method can be made more efficient (because removal of each
-  // inequality leads to much shifting/copying in the underlying buffer).
-  for (auto nbIndex : llvm::reverse(nbIneqIndices))
-    removeInequality(nbIndex);
-  for (auto nbIndex : llvm::reverse(nbEqIndices))
-    removeEquality(nbIndex);
-}
-
 AffineMap mlir::alignAffineMapWithValues(AffineMap map, ValueRange operands,
                                          ValueRange dims, ValueRange syms,
                                          SmallVector<Value> *newSyms) {
   assert(operands.size() == map.getNumInputs() &&
          "expected same number of operands and map inputs");
   MLIRContext *ctx = map.getContext();
   Builder builder(ctx);
   SmallVector<AffineExpr> dimReplacements(map.getNumDims(), {});

mlir/lib/Analysis/CMakeLists.txt

 set(LLVM_OPTIONAL_SOURCES
   AliasAnalysis.cpp
   AffineAnalysis.cpp
   AffineStructures.cpp
   BufferViewFlowAnalysis.cpp
   CallGraph.cpp
   DataFlowAnalysis.cpp
   DataLayoutAnalysis.cpp
-  LinearTransform.cpp
   Liveness.cpp
   LoopAnalysis.cpp
   NestedMatcher.cpp
   NumberOfExecutions.cpp
   PresburgerSet.cpp
   SliceAnalysis.cpp
   Utils.cpp
 
add_mlir_library(MLIRAnalysis
   MLIRInferTypeOpInterface
   MLIRSideEffectInterfaces
   MLIRViewLikeInterface
   )
 
 add_mlir_library(MLIRLoopAnalysis
   AffineAnalysis.cpp
   AffineStructures.cpp
-  LinearTransform.cpp
   LoopAnalysis.cpp
   NestedMatcher.cpp
   PresburgerSet.cpp
   Utils.cpp
 
   ADDITIONAL_HEADER_DIRS
   ${MLIR_MAIN_INCLUDE_DIR}/mlir/Analysis
 

mlir/lib/Analysis/Presburger/CMakeLists.txt

 add_mlir_library(MLIRPresburger
   IntegerPolyhedron.cpp
+  LinearTransform.cpp
   Matrix.cpp
   Simplex.cpp
   Utils.cpp
 
   DEPENDS
   MLIRBuiltinLocationAttributesIncGen
 
   LINK_LIBS PUBLIC
   MLIRIR
   MLIRSupport
   )

mlir/lib/Analysis/Presburger/IntegerPolyhedron.cpp

 //===- IntegerPolyhedron.cpp - MLIR IntegerPolyhedron 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
 //
 //===----------------------------------------------------------------------===//
 //
 // A class to represent an integer polyhedron.
 //
 //===----------------------------------------------------------------------===//
 
 #include "mlir/Analysis/Presburger/IntegerPolyhedron.h"
+#include "mlir/Analysis/Presburger/LinearTransform.h"
+#include "mlir/Analysis/Presburger/Simplex.h"
+#include "mlir/Analysis/Presburger/Utils.h"
+#include "mlir/Support/LLVM.h"
+#include "mlir/Support/MathExtras.h"
+#include "llvm/ADT/STLExtras.h"
+#include "llvm/Support/Debug.h"
+
+#define DEBUG_TYPE "presburger"
 
 using namespace mlir;
+using llvm::SmallDenseMap;
+using llvm::SmallDenseSet;
 
 std::unique_ptr<IntegerPolyhedron> IntegerPolyhedron::clone() const {
   return std::make_unique<IntegerPolyhedron>(*this);
 }
 
 void IntegerPolyhedron::reset(unsigned numReservedInequalities,
                               unsigned numReservedEqualities,
                               unsigned newNumReservedCols, unsigned newNumDims,
void IntegerPolyhedron::getLowerAndUpperBoundIndices(
   for (unsigned r = 0, e = getNumEqualities(); r < e; r++) {
     if (atEq(r, pos) == 0)
       continue;
     if (containsConstraintDependentOnRange(r, /*isEq=*/true))
       continue;
     eqIndices->push_back(r);
   }
 }
+
+// Searches for a constraint with a non-zero coefficient at `colIdx` in
+// equality (isEq=true) or inequality (isEq=false) constraints.
+// Returns true and sets row found in search in `rowIdx`, false otherwise.
+bool IntegerPolyhedron::findConstraintWithNonZeroAt(unsigned colIdx, bool isEq,
+                                                    unsigned *rowIdx) const {
+  assert(colIdx < getNumCols() && "position out of bounds");
+  auto at = [&](unsigned rowIdx) -> int64_t {
+    return isEq ? atEq(rowIdx, colIdx) : atIneq(rowIdx, colIdx);
+  };
+  unsigned e = isEq ? getNumEqualities() : getNumInequalities();
+  for (*rowIdx = 0; *rowIdx < e; ++(*rowIdx)) {
+    if (at(*rowIdx) != 0) {
+      return true;
+    }
+  }
+  return false;
+}
+
+// Normalizes the coefficient values across all columns in `rowIdx` by their
+// GCD in equality or inequality constraints as specified by `isEq`.
+template <bool isEq>
+static void normalizeConstraintByGCD(IntegerPolyhedron *constraints,
+                                     unsigned rowIdx) {
+  auto at = [&](unsigned colIdx) -> int64_t {
+    return isEq ? constraints->atEq(rowIdx, colIdx)
+                : constraints->atIneq(rowIdx, colIdx);
+  };
+  uint64_t gcd = std::abs(at(0));
+  for (unsigned j = 1, e = constraints->getNumCols(); j < e; ++j) {
+    gcd = llvm::GreatestCommonDivisor64(gcd, std::abs(at(j)));
+  }
+  if (gcd > 0 && gcd != 1) {
+    for (unsigned j = 0, e = constraints->getNumCols(); j < e; ++j) {
+      int64_t v = at(j) / static_cast<int64_t>(gcd);
+      isEq ? constraints->atEq(rowIdx, j) = v
+           : constraints->atIneq(rowIdx, j) = v;
+    }
+  }
+}
+
+void IntegerPolyhedron::normalizeConstraintsByGCD() {
+  for (unsigned i = 0, e = getNumEqualities(); i < e; ++i) {
+    normalizeConstraintByGCD</*isEq=*/true>(this, i);
+  }
+  for (unsigned i = 0, e = getNumInequalities(); i < e; ++i) {
+    normalizeConstraintByGCD</*isEq=*/false>(this, i);
+  }
+}
+
+bool IntegerPolyhedron::hasConsistentState() const {
+  if (!inequalities.hasConsistentState())
+    return false;
+  if (!equalities.hasConsistentState())
+    return false;
+
+  // Catches errors where numDims, numSymbols, numIds aren't consistent.
+  if (numDims > numIds || numSymbols > numIds || numDims + numSymbols > numIds)
+    return false;
+
+  return true;
+}
+
+bool IntegerPolyhedron::hasInvalidConstraint() const {
+  assert(hasConsistentState());
+  auto check = [&](bool isEq) -> bool {
+    unsigned numCols = getNumCols();
+    unsigned numRows = isEq ? getNumEqualities() : getNumInequalities();
+    for (unsigned i = 0, e = numRows; i < e; ++i) {
+      unsigned j;
+      for (j = 0; j < numCols - 1; ++j) {
+        int64_t v = isEq ? atEq(i, j) : atIneq(i, j);
+        // Skip rows with non-zero variable coefficients.
+        if (v != 0)
+          break;
+      }
+      if (j < numCols - 1) {
+        continue;
+      }
+      // Check validity of constant term at 'numCols - 1' w.r.t 'isEq'.
+      // Example invalid constraints include: '1 == 0' or '-1 >= 0'
+      int64_t v = isEq ? atEq(i, numCols - 1) : atIneq(i, numCols - 1);
+      if ((isEq && v != 0) || (!isEq && v < 0)) {
+        return true;
+      }
+    }
+    return false;
+  };
+  if (check(/*isEq=*/true))
+    return true;
+  return check(/*isEq=*/false);
+}
+
+/// Eliminate identifier from constraint at `rowIdx` based on coefficient at
+/// pivotRow, pivotCol. Columns in range [elimColStart, pivotCol) will not be
+/// updated as they have already been eliminated.
+static void eliminateFromConstraint(IntegerPolyhedron *constraints,
+                                    unsigned rowIdx, unsigned pivotRow,
+                                    unsigned pivotCol, unsigned elimColStart,
+                                    bool isEq) {
+  // Skip if equality 'rowIdx' if same as 'pivotRow'.
+  if (isEq && rowIdx == pivotRow)
+    return;
+  auto at = [&](unsigned i, unsigned j) -> int64_t {
+    return isEq ? constraints->atEq(i, j) : constraints->atIneq(i, j);
+  };
+  int64_t leadCoeff = at(rowIdx, pivotCol);
+  // Skip if leading coefficient at 'rowIdx' is already zero.
+  if (leadCoeff == 0)
+    return;
+  int64_t pivotCoeff = constraints->atEq(pivotRow, pivotCol);
+  int64_t sign = (leadCoeff * pivotCoeff > 0) ? -1 : 1;
+  int64_t lcm = mlir::lcm(pivotCoeff, leadCoeff);
+  int64_t pivotMultiplier = sign * (lcm / std::abs(pivotCoeff));
+  int64_t rowMultiplier = lcm / std::abs(leadCoeff);
+
+  unsigned numCols = constraints->getNumCols();
+  for (unsigned j = 0; j < numCols; ++j) {
+    // Skip updating column 'j' if it was just eliminated.
+    if (j >= elimColStart && j < pivotCol)
+      continue;
+    int64_t v = pivotMultiplier * constraints->atEq(pivotRow, j) +
+                rowMultiplier * at(rowIdx, j);
+    isEq ? constraints->atEq(rowIdx, j) = v
+         : constraints->atIneq(rowIdx, j) = v;
+  }
+}
+
+/// Returns the position of the identifier that has the minimum <number of lower
+/// bounds> times <number of upper bounds> from the specified range of
+/// identifiers [start, end). It is often best to eliminate in the increasing
+/// order of these counts when doing Fourier-Motzkin elimination since FM adds
+/// that many new constraints.
+static unsigned getBestIdToEliminate(const IntegerPolyhedron &cst,
+                                     unsigned start, unsigned end) {
+  assert(start < cst.getNumIds() && end < cst.getNumIds() + 1);
+
+  auto getProductOfNumLowerUpperBounds = [&](unsigned pos) {
+    unsigned numLb = 0;
+    unsigned numUb = 0;
+    for (unsigned r = 0, e = cst.getNumInequalities(); r < e; r++) {
+      if (cst.atIneq(r, pos) > 0) {
+        ++numLb;
+      } else if (cst.atIneq(r, pos) < 0) {
+        ++numUb;
+      }
+    }
+    return numLb * numUb;
+  };
+
+  unsigned minLoc = start;
+  unsigned min = getProductOfNumLowerUpperBounds(start);
+  for (unsigned c = start + 1; c < end; c++) {
+    unsigned numLbUbProduct = getProductOfNumLowerUpperBounds(c);
+    if (numLbUbProduct < min) {
+      min = numLbUbProduct;
+      minLoc = c;
+    }
+  }
+  return minLoc;
+}
+
+// Checks for emptiness of the set by eliminating identifiers successively and
+// using the GCD test (on all equality constraints) and checking for trivially
+// invalid constraints. Returns 'true' if the constraint system is found to be
+// empty; false otherwise.
+bool IntegerPolyhedron::isEmpty() const {
+  if (isEmptyByGCDTest() || hasInvalidConstraint())
+    return true;
+
+  IntegerPolyhedron tmpCst(*this);
+
+  // First, eliminate as many local variables as possible using equalities.
+  tmpCst.removeRedundantLocalVars();
+  if (tmpCst.isEmptyByGCDTest() || tmpCst.hasInvalidConstraint())
+    return true;
+
+  // Eliminate as many identifiers as possible using Gaussian elimination.
+  unsigned currentPos = 0;
+  while (currentPos < tmpCst.getNumIds()) {
+    tmpCst.gaussianEliminateIds(currentPos, tmpCst.getNumIds());
+    ++currentPos;
+    // We check emptiness through trivial checks after eliminating each ID to
+    // detect emptiness early. Since the checks isEmptyByGCDTest() and
+    // hasInvalidConstraint() are linear time and single sweep on the constraint
+    // buffer, this appears reasonable - but can optimize in the future.
+    if (tmpCst.hasInvalidConstraint() || tmpCst.isEmptyByGCDTest())
+      return true;
+  }
+
+  // Eliminate the remaining using FM.
+  for (unsigned i = 0, e = tmpCst.getNumIds(); i < e; i++) {
+    tmpCst.fourierMotzkinEliminate(
+        getBestIdToEliminate(tmpCst, 0, tmpCst.getNumIds()));
+    // Check for a constraint explosion. This rarely happens in practice, but
+    // this check exists as a safeguard against improperly constructed
+    // constraint systems or artificially created arbitrarily complex systems
+    // that aren't the intended use case for IntegerPolyhedron. This is
+    // needed since FM has a worst case exponential complexity in theory.
+    if (tmpCst.getNumConstraints() >= kExplosionFactor * getNumIds()) {
+      LLVM_DEBUG(llvm::dbgs() << "FM constraint explosion detected\n");
+      return false;
+    }
+
+    // FM wouldn't have modified the equalities in any way. So no need to again
+    // run GCD test. Check for trivial invalid constraints.
+    if (tmpCst.hasInvalidConstraint())
+      return true;
+  }
+  return false;
+}
+
+// Runs the GCD test on all equality constraints. Returns 'true' if this test
+// fails on any equality. Returns 'false' otherwise.
+// This test can be used to disprove the existence of a solution. If it returns
+// true, no integer solution to the equality constraints can exist.
+//
+// GCD test definition:
+//
+// The equality constraint:
+//
+//  c_1*x_1 + c_2*x_2 + ... + c_n*x_n = c_0
+//
+// has an integer solution iff:
+//
+//  GCD of c_1, c_2, ..., c_n divides c_0.
+//
+bool IntegerPolyhedron::isEmptyByGCDTest() const {
+  assert(hasConsistentState());
+  unsigned numCols = getNumCols();
+  for (unsigned i = 0, e = getNumEqualities(); i < e; ++i) {
+    uint64_t gcd = std::abs(atEq(i, 0));
+    for (unsigned j = 1; j < numCols - 1; ++j) {
+      gcd = llvm::GreatestCommonDivisor64(gcd, std::abs(atEq(i, j)));
+    }
+    int64_t v = std::abs(atEq(i, numCols - 1));
+    if (gcd > 0 && (v % gcd != 0)) {
+      return true;
+    }
+  }
+  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.
+//
+// It is sufficient to check the perpendiculars of the constraints, as the set
+// of perpendiculars which are bounded must span all bounded directions.
+Matrix IntegerPolyhedron::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.
+  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 IntegerPolyhedron &fac, unsigned eqIndex,
+                          unsigned numDims) {
+  for (unsigned e = fac.getNumIds(), j = e - numDims; j < e; ++j)
+    if (fac.atEq(eqIndex, j) != 0)
+      return true;
+  return false;
+}
+bool ineqInvolvesSuffixDims(const IntegerPolyhedron &fac, unsigned ineqIndex,
+                            unsigned numDims) {
+  for (unsigned e = fac.getNumIds(), j = e - numDims; j < e; ++j)
+    if (fac.atIneq(ineqIndex, j) != 0)
+      return true;
+  return false;
+}
+
+void removeConstraintsInvolvingSuffixDims(IntegerPolyhedron &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);
+}
+
+bool IntegerPolyhedron::isIntegerEmpty() const {
+  return !findIntegerSample().hasValue();
+}
+
+/// 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 infinity 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 B by removing all constraints that involve
+/// the unbounded dimensions and then deleting the unbounded dimensions. Note
+/// that B is a Bounded set.
+///
+/// 3) Try to obtain a sample from B using the GBR sampling
+/// algorithm. If no sample is found, return that S is empty.
+///
+/// 4) Otherwise, substitute the obtained sample into S*T to obtain a set
+/// C. C is a full-dimensional Cone and always contains a sample.
+///
+/// 5) Obtain an integer sample from C.
+///
+/// 6) Return T*v, where v is the concatenation of the samples from B and C.
+///
+/// The following is a sketch of a proof that
+/// a) If the algorithm returns empty, then S is empty.
+/// b) If the algorithm returns a sample, it is a valid sample in S.
+///
+/// The algorithm returns empty only if B is empty, in which case S*T is
+/// certainly empty since B was obtained by removing constraints and then
+/// deleting unconstrained dimensions from S*T. Since T is unimodular, a vector
+/// v is in S*T iff T*v is in S. So in this case, since
+/// S*T is empty, S is empty too.
+///
+/// Otherwise, the algorithm substitutes the sample from B into S*T. All the
+/// constraints of S*T that did not involve unbounded dimensions are satisfied
+/// by this substitution. All dimensions in the linear span of the dimensions
+/// outside the prefix are unbounded in S*T (step 1). Substituting values for
+/// the bounded dimensions cannot make these dimensions bounded, and these are
+/// the only remaining dimensions in C, so C is unbounded along every vector (in
+/// the positive or negative direction, or both). C is hence a full-dimensional
+/// cone and therefore always contains an integer point.
+///
+/// Concatenating the samples from B and C gives a sample v in S*T, so the
+/// returned sample T*v is a sample in S.
+Optional<SmallVector<int64_t, 8>> IntegerPolyhedron::findIntegerSample() const {
+  // First, try the GCD test heuristic.
+  if (isEmptyByGCDTest())
+    return {};
+
+  Simplex simplex(*this);
+  if (simplex.isEmpty())
+    return {};
+
+  // For a bounded set, we directly call into the GBR sampling algorithm.
+  if (!simplex.isUnbounded())
+    return simplex.findIntegerSample();
+
+  // 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<unsigned, LinearTransform> result =
+      LinearTransform::makeTransformToColumnEchelon(std::move(m));
+  const LinearTransform &transform = result.second;
+  // 1) Apply T to S to obtain S*T.
+  IntegerPolyhedron transformedSet = transform.applyTo(*this);
+
+  // 2) Remove the unbounded dimensions and constraints involving them to
+  // obtain a bounded set.
+  IntegerPolyhedron boundedSet = transformedSet;
+  unsigned numBoundedDims = result.first;
+  unsigned numUnboundedDims = getNumIds() - numBoundedDims;
+  removeConstraintsInvolvingSuffixDims(boundedSet, numUnboundedDims);
+  boundedSet.removeIdRange(numBoundedDims, boundedSet.getNumIds());
+
+  // 3) Try to obtain a sample from the bounded set.
+  Optional<SmallVector<int64_t, 8>> boundedSample =
+      Simplex(boundedSet).findIntegerSample();
+  if (!boundedSample)
+    return {};
+  assert(boundedSet.containsPoint(*boundedSample) &&
+         "Simplex returned an invalid sample!");
+
+  // 4) Substitute the values of the bounded dimensions into S*T to obtain a
+  // full-dimensional cone, which necessarily contains an integer sample.
+  transformedSet.setAndEliminate(0, *boundedSample);
+  IntegerPolyhedron &cone = transformedSet;
+
+  // 5) Obtain an integer sample from the cone.
+  //
+  // We shrink the cone such that for any rational point in the shrunken cone,
+  // rounding up each of the point's coordinates produces a point that still
+  // lies in the original cone.
+  //
+  // Rounding up a point x adds a number e_i in [0, 1) to each coordinate x_i.
+  // For each inequality sum_i a_i x_i + c >= 0 in the original cone, the
+  // shrunken cone will have the inequality tightened by some amount s, such
+  // that if x satisfies the shrunken cone's tightened inequality, then x + e
+  // satisfies the original inequality, i.e.,
+  //
+  // sum_i a_i x_i + c + s >= 0 implies sum_i a_i (x_i + e_i) + c >= 0
+  //
+  // for any e_i values in [0, 1). In fact, we will handle the slightly more
+  // general case where e_i can be in [0, 1]. For example, consider the
+  // inequality 2x_1 - 3x_2 - 7x_3 - 6 >= 0, and let x = (3, 0, 0). How low
+  // could the LHS go if we added a number in [0, 1] to each coordinate? The LHS
+  // is minimized when we add 1 to the x_i with negative coefficient a_i and
+  // keep the other x_i the same. In the example, we would get x = (3, 1, 1),
+  // changing the value of the LHS by -3 + -7 = -10.
+  //
+  // In general, the value of the LHS can change by at most the sum of the
+  // negative a_i, so we accomodate this by shifting the inequality by this
+  // amount for the shrunken cone.
+  for (unsigned i = 0, e = cone.getNumInequalities(); i < e; ++i) {
+    for (unsigned j = 0; j < cone.numIds; ++j) {
+      int64_t coeff = cone.atIneq(i, j);
+      if (coeff < 0)
+        cone.atIneq(i, cone.numIds) += coeff;
+    }
+  }
+
+  // Obtain an integer sample in the cone by rounding up a rational point from
+  // the shrunken cone. Shrinking the cone amounts to shifting its apex
+  // "inwards" without changing its "shape"; the shrunken cone is still a
+  // full-dimensional cone and is hence non-empty.
+  Simplex shrunkenConeSimplex(cone);
+  assert(!shrunkenConeSimplex.isEmpty() && "Shrunken cone cannot be empty!");
+  SmallVector<Fraction, 8> shrunkenConeSample =
+      shrunkenConeSimplex.getRationalSample();
+
+  SmallVector<int64_t, 8> coneSample(llvm::map_range(shrunkenConeSample, ceil));
+
+  // 6) Return transform * concat(boundedSample, coneSample).
+  SmallVector<int64_t, 8> &sample = boundedSample.getValue();
+  sample.append(coneSample.begin(), coneSample.end());
+  return transform.preMultiplyColumn(sample);
+}
+
+/// Helper to evaluate an affine expression at a point.
+/// The expression is a list of coefficients for the dimensions followed by the
+/// constant term.
+static int64_t valueAt(ArrayRef<int64_t> expr, ArrayRef<int64_t> point) {
+  assert(expr.size() == 1 + point.size() &&
+         "Dimensionalities of point and expression don't match!");
+  int64_t value = expr.back();
+  for (unsigned i = 0; i < point.size(); ++i)
+    value += expr[i] * point[i];
+  return value;
+}
+
+/// A point satisfies an equality iff the value of the equality at the
+/// expression is zero, and it satisfies an inequality iff the value of the
+/// inequality at that point is non-negative.
+bool IntegerPolyhedron::containsPoint(ArrayRef<int64_t> point) const {
+  for (unsigned i = 0, e = getNumEqualities(); i < e; ++i) {
+    if (valueAt(getEquality(i), point) != 0)
+      return false;
+  }
+  for (unsigned i = 0, e = getNumInequalities(); i < e; ++i) {
+    if (valueAt(getInequality(i), point) < 0)
+      return false;
+  }
+  return true;
+}
+
+void IntegerPolyhedron::getLocalReprs(
+    std::vector<llvm::Optional<std::pair<unsigned, unsigned>>> &repr) const {
+  std::vector<SmallVector<int64_t, 8>> dividends(getNumLocalIds());
+  SmallVector<unsigned, 4> denominators(getNumLocalIds());
+  getLocalReprs(dividends, denominators, repr);
+}
+
+void IntegerPolyhedron::getLocalReprs(
+    std::vector<SmallVector<int64_t, 8>> &dividends,
+    SmallVector<unsigned, 4> &denominators) const {
+  std::vector<llvm::Optional<std::pair<unsigned, unsigned>>> repr(
+      getNumLocalIds());
+  getLocalReprs(dividends, denominators, repr);
+}
+
+void IntegerPolyhedron::getLocalReprs(
+    std::vector<SmallVector<int64_t, 8>> &dividends,
+    SmallVector<unsigned, 4> &denominators,
+    std::vector<llvm::Optional<std::pair<unsigned, unsigned>>> &repr) const {
+
+  repr.resize(getNumLocalIds());
+  dividends.resize(getNumLocalIds());
+  denominators.resize(getNumLocalIds());
+
+  SmallVector<bool, 8> foundRepr(getNumIds(), false);
+  for (unsigned i = 0, e = getNumDimAndSymbolIds(); i < e; ++i)
+    foundRepr[i] = true;
+
+  unsigned divOffset = getNumDimAndSymbolIds();
+  bool changed;
+  do {
+    // Each time changed is true, at end of this iteration, one or more local
+    // vars have been detected as floor divs.
+    changed = false;
+    for (unsigned i = 0, e = getNumLocalIds(); i < e; ++i) {
+      if (!foundRepr[i + divOffset]) {
+        if (auto res = presburger_utils::computeSingleVarRepr(
+                *this, foundRepr, divOffset + i, dividends[i],
+                denominators[i])) {
+          foundRepr[i + divOffset] = true;
+          repr[i] = res;
+          changed = true;
+        }
+      }
+    }
+  } while (changed);
+
+  // Set 0 denominator for identifiers for which no division representation
+  // could be found.
+  for (unsigned i = 0, e = repr.size(); i < e; ++i)
+    if (!repr[i].hasValue())
+      denominators[i] = 0;
+}
+
+/// Tightens inequalities given that we are dealing with integer spaces. This is
+/// analogous to the GCD test but applied to inequalities. The constant term can
+/// be reduced to the preceding multiple of the GCD of the coefficients, i.e.,
+///  64*i - 100 >= 0  =>  64*i - 128 >= 0 (since 'i' is an integer). This is a
+/// fast method - linear in the number of coefficients.
+// Example on how this affects practical cases: consider the scenario:
+// 64*i >= 100, j = 64*i; without a tightening, elimination of i would yield
+// j >= 100 instead of the tighter (exact) j >= 128.
+void IntegerPolyhedron::gcdTightenInequalities() {
+  unsigned numCols = getNumCols();
+  for (unsigned i = 0, e = getNumInequalities(); i < e; ++i) {
+    uint64_t gcd = std::abs(atIneq(i, 0));
+    for (unsigned j = 1; j < numCols - 1; ++j) {
+      gcd = llvm::GreatestCommonDivisor64(gcd, std::abs(atIneq(i, j)));
+    }
+    if (gcd > 0 && gcd != 1) {
+      int64_t gcdI = static_cast<int64_t>(gcd);
+      // Tighten the constant term and normalize the constraint by the GCD.
+      atIneq(i, numCols - 1) = mlir::floorDiv(atIneq(i, numCols - 1), gcdI);
+      for (unsigned j = 0, e = numCols - 1; j < e; ++j)
+        atIneq(i, j) /= gcdI;
+    }
+  }
+}
+
+// Eliminates all identifier variables in column range [posStart, posLimit).
+// Returns the number of variables eliminated.
+unsigned IntegerPolyhedron::gaussianEliminateIds(unsigned posStart,
+                                                 unsigned posLimit) {
+  // Return if identifier positions to eliminate are out of range.
+  assert(posLimit <= numIds);
+  assert(hasConsistentState());
+
+  if (posStart >= posLimit)
+    return 0;
+
+  gcdTightenInequalities();
+
+  unsigned pivotCol = 0;
+  for (pivotCol = posStart; pivotCol < posLimit; ++pivotCol) {
+    // Find a row which has a non-zero coefficient in column 'j'.
+    unsigned pivotRow;
+    if (!findConstraintWithNonZeroAt(pivotCol, /*isEq=*/true, &pivotRow)) {
+      // No pivot row in equalities with non-zero at 'pivotCol'.
+      if (!findConstraintWithNonZeroAt(pivotCol, /*isEq=*/false, &pivotRow)) {
+        // If inequalities are also non-zero in 'pivotCol', it can be
+        // eliminated.
+        continue;
+      }
+      break;
+    }
+
+    // Eliminate identifier at 'pivotCol' from each equality row.
+    for (unsigned i = 0, e = getNumEqualities(); i < e; ++i) {
+      eliminateFromConstraint(this, i, pivotRow, pivotCol, posStart,
+                              /*isEq=*/true);
+      normalizeConstraintByGCD</*isEq=*/true>(this, i);
+    }
+
+    // Eliminate identifier at 'pivotCol' from each inequality row.
+    for (unsigned i = 0, e = getNumInequalities(); i < e; ++i) {
+      eliminateFromConstraint(this, i, pivotRow, pivotCol, posStart,
+                              /*isEq=*/false);
+      normalizeConstraintByGCD</*isEq=*/false>(this, i);
+    }
+    removeEquality(pivotRow);
+    gcdTightenInequalities();
+  }
+  // Update position limit based on number eliminated.
+  posLimit = pivotCol;
+  // Remove eliminated columns from all constraints.
+  removeIdRange(posStart, posLimit);
+  return posLimit - posStart;
+}
+
+// Fills an inequality row with the value 'val'.
+static inline void fillInequality(IntegerPolyhedron *cst, unsigned r,
+                                  int64_t val) {
+  for (unsigned c = 0, f = cst->getNumCols(); c < f; c++) {
+    cst->atIneq(r, c) = val;
+  }
+}
+
+// Negates an inequality.
+static inline void negateInequality(IntegerPolyhedron *cst, unsigned r) {
+  for (unsigned c = 0, f = cst->getNumCols(); c < f; c++) {
+    cst->atIneq(r, c) = -cst->atIneq(r, c);
+  }
+}
+
+// A more complex check to eliminate redundant inequalities. Uses FourierMotzkin
+// to check if a constraint is redundant.
+void IntegerPolyhedron::removeRedundantInequalities() {
+  SmallVector<bool, 32> redun(getNumInequalities(), false);
+  // To check if an inequality is redundant, we replace the inequality by its
+  // complement (for eg., i - 1 >= 0 by i <= 0), and check if the resulting
+  // system is empty. If it is, the inequality is redundant.
+  IntegerPolyhedron tmpCst(*this);
+  for (unsigned r = 0, e = getNumInequalities(); r < e; r++) {
+    // Change the inequality to its complement.
+    negateInequality(&tmpCst, r);
+    tmpCst.atIneq(r, tmpCst.getNumCols() - 1)--;
+    if (tmpCst.isEmpty()) {
+      redun[r] = true;
+      // Zero fill the redundant inequality.
+      fillInequality(this, r, /*val=*/0);
+      fillInequality(&tmpCst, r, /*val=*/0);
+    } else {
+      // Reverse the change (to avoid recreating tmpCst each time).
+      tmpCst.atIneq(r, tmpCst.getNumCols() - 1)++;
+      negateInequality(&tmpCst, r);
+    }
+  }
+
+  // Scan to get rid of all rows marked redundant, in-place.
+  auto copyRow = [&](unsigned src, unsigned dest) {
+    if (src == dest)
+      return;
+    for (unsigned c = 0, e = getNumCols(); c < e; c++) {
+      atIneq(dest, c) = atIneq(src, c);
+    }
+  };
+  unsigned pos = 0;
+  for (unsigned r = 0, e = getNumInequalities(); r < e; r++) {
+    if (!redun[r])
+      copyRow(r, pos++);
+  }
+  inequalities.resizeVertically(pos);
+}
+
+// A more complex check to eliminate redundant inequalities and equalities. Uses
+// Simplex to check if a constraint is redundant.
+void IntegerPolyhedron::removeRedundantConstraints() {
+  // First, we run gcdTightenInequalities. This allows us to catch some
+  // constraints which are not redundant when considering rational solutions
+  // but are redundant in terms of integer solutions.
+  gcdTightenInequalities();
+  Simplex simplex(*this);
+  simplex.detectRedundant();
+
+  auto copyInequality = [&](unsigned src, unsigned dest) {
+    if (src == dest)
+      return;
+    for (unsigned c = 0, e = getNumCols(); c < e; c++)
+      atIneq(dest, c) = atIneq(src, c);
+  };
+  unsigned pos = 0;
+  unsigned numIneqs = getNumInequalities();
+  // Scan to get rid of all inequalities marked redundant, in-place. In Simplex,
+  // the first constraints added are the inequalities.
+  for (unsigned r = 0; r < numIneqs; r++) {
+    if (!simplex.isMarkedRedundant(r))
+      copyInequality(r, pos++);
+  }
+  inequalities.resizeVertically(pos);
+
+  // Scan to get rid of all equalities marked redundant, in-place. In Simplex,
+  // after the inequalities, a pair of constraints for each equality is added.
+  // An equality is redundant if both the inequalities in its pair are
+  // redundant.
+  auto copyEquality = [&](unsigned src, unsigned dest) {
+    if (src == dest)
+      return;
+    for (unsigned c = 0, e = getNumCols(); c < e; c++)
+      atEq(dest, c) = atEq(src, c);
+  };
+  pos = 0;
+  for (unsigned r = 0, e = getNumEqualities(); r < e; r++) {
+    if (!(simplex.isMarkedRedundant(numIneqs + 2 * r) &&
+          simplex.isMarkedRedundant(numIneqs + 2 * r + 1)))
+      copyEquality(r, pos++);
+  }
+  equalities.resizeVertically(pos);
+}
+
+/// Eliminate `pos2^th` local identifier, replacing its every instance with
+/// `pos1^th` local identifier. This function is intended to be used to remove
+/// redundancy when local variables at position `pos1` and `pos2` are restricted
+/// to have the same value.
+static void eliminateRedundantLocalId(IntegerPolyhedron &fac, unsigned pos1,
+                                      unsigned pos2) {
+
+  assert(pos1 < fac.getNumLocalIds() && "Invalid local id position");
+  assert(pos2 < fac.getNumLocalIds() && "Invalid local id position");
+
+  unsigned localOffset = fac.getNumDimAndSymbolIds();
+  pos1 += localOffset;
+  pos2 += localOffset;
+  for (unsigned i = 0, e = fac.getNumInequalities(); i < e; ++i)
+    fac.atIneq(i, pos1) += fac.atIneq(i, pos2);
+  for (unsigned i = 0, e = fac.getNumEqualities(); i < e; ++i)
+    fac.atEq(i, pos1) += fac.atEq(i, pos2);
+  fac.removeId(pos2);
+}
+
+/// Adds additional local ids to the sets such that they both have the union
+/// of the local ids in each set, without changing the set of points that
+/// lie in `this` and `other`.
+///
+/// To detect local ids that always take the same in both sets, each local id is
+/// represented as a floordiv with constant denominator in terms of other ids.
+/// After extracting these divisions, local ids with the same division
+/// representation are considered duplicate and are merged. It is possible that
+/// division representation for some local id cannot be obtained, and thus these
+/// local ids are not considered for detecting duplicates.
+void IntegerPolyhedron::mergeLocalIds(IntegerPolyhedron &other) {
+  assert(getNumDimIds() == other.getNumDimIds() &&
+         "Number of dimension ids should match");
+  assert(getNumSymbolIds() == other.getNumSymbolIds() &&
+         "Number of symbol ids should match");
+
+  IntegerPolyhedron &facA = *this;
+  IntegerPolyhedron &facB = other;
+
+  // Merge local ids of facA and facB without using division information,
+  // i.e. append local ids of `facB` to `facA` and insert local ids of `facA`
+  // to `facB` at start of its local ids.
+  unsigned initLocals = facA.getNumLocalIds();
+  insertLocalId(facA.getNumLocalIds(), facB.getNumLocalIds());
+  facB.insertLocalId(0, initLocals);
+
+  // Get division representations from each FAC.
+  std::vector<SmallVector<int64_t, 8>> divsA, divsB;
+  SmallVector<unsigned, 4> denomsA, denomsB;
+  facA.getLocalReprs(divsA, denomsA);
+  facB.getLocalReprs(divsB, denomsB);
+
+  // Copy division information for facB into `divsA` and `denomsA`, so that
+  // these have the combined division information of both FACs. Since newly
+  // added local variables in facA and facB have no constraints, they will not
+  // have any division representation.
+  std::copy(divsB.begin() + initLocals, divsB.end(),
+            divsA.begin() + initLocals);
+  std::copy(denomsB.begin() + initLocals, denomsB.end(),
+            denomsA.begin() + initLocals);
+
+  // Find and merge duplicate divisions.
+  // TODO: Add division normalization to support divisions that differ by
+  // a constant.
+  // TODO: Add division ordering such that a division representation for local
+  // identifier at position `i` only depends on local identifiers at position <
+  // `i`. This would make sure that all divisions depending on other local
+  // variables that can be merged, are merged.
+  unsigned localOffset = getIdKindOffset(IdKind::Local);
+  for (unsigned i = 0; i < divsA.size(); ++i) {
+    // Check if a division representation exists for the `i^th` local id.
+    if (denomsA[i] == 0)
+      continue;
+    // Check if a division exists which is a duplicate of the division at `i`.
+    for (unsigned j = i + 1; j < divsA.size(); ++j) {
+      // Check if a division representation exists for the `j^th` local id.
+      if (denomsA[j] == 0)
+        continue;
+      // Check if the denominators match.
+      if (denomsA[i] != denomsA[j])
+        continue;
+      // Check if the representations are equal.
+      if (divsA[i] != divsA[j])
+        continue;
+
+      // Merge divisions at position `j` into division at position `i`.
+      eliminateRedundantLocalId(facA, i, j);
+      eliminateRedundantLocalId(facB, i, j);
+      for (unsigned k = 0, g = divsA.size(); k < g; ++k) {
+        SmallVector<int64_t, 8> &div = divsA[k];
+        if (denomsA[k] != 0) {
+          div[localOffset + i] += div[localOffset + j];
+          div.erase(div.begin() + localOffset + j);
+        }
+      }
+
+      divsA.erase(divsA.begin() + j);
+      denomsA.erase(denomsA.begin() + j);
+      // Since `j` can never be zero, we do not need to worry about overflows.
+      --j;
+    }
+  }
+}
+
+/// Removes local variables using equalities. Each equality is checked if it
+/// can be reduced to the form: `e = affine-expr`, where `e` is a local
+/// variable and `affine-expr` is an affine expression not containing `e`.
+/// If an equality satisfies this form, the local variable is replaced in
+/// each constraint and then removed. The equality used to replace this local
+/// variable is also removed.
+void IntegerPolyhedron::removeRedundantLocalVars() {
+  // Normalize the equality constraints to reduce coefficients of local
+  // variables to 1 wherever possible.
+  for (unsigned i = 0, e = getNumEqualities(); i < e; ++i)
+    normalizeConstraintByGCD</*isEq=*/true>(this, i);
+
+  while (true) {
+    unsigned i, e, j, f;
+    for (i = 0, e = getNumEqualities(); i < e; ++i) {
+      // Find a local variable to eliminate using ith equality.
+      for (j = getNumDimAndSymbolIds(), f = getNumIds(); j < f; ++j)
+        if (std::abs(atEq(i, j)) == 1)
+          break;
+
+      // Local variable can be eliminated using ith equality.
+      if (j < f)
+        break;
+    }
+
+    // No equality can be used to eliminate a local variable.
+    if (i == e)
+      break;
+
+    // Use the ith equality to simplify other equalities. If any changes
+    // are made to an equality constraint, it is normalized by GCD.
+    for (unsigned k = 0, t = getNumEqualities(); k < t; ++k) {
+      if (atEq(k, j) != 0) {
+        eliminateFromConstraint(this, k, i, j, j, /*isEq=*/true);
+        normalizeConstraintByGCD</*isEq=*/true>(this, k);
+      }
+    }
+
+    // Use the ith equality to simplify inequalities.
+    for (unsigned k = 0, t = getNumInequalities(); k < t; ++k)
+      eliminateFromConstraint(this, k, i, j, j, /*isEq=*/false);
+
+    // Remove the ith equality and the found local variable.
+    removeId(j);
+    removeEquality(i);
+  }
+}
+
+void IntegerPolyhedron::convertDimToLocal(unsigned dimStart,
+                                          unsigned dimLimit) {
+  assert(dimLimit <= getNumDimIds() && "Invalid dim pos range");
+
+  if (dimStart >= dimLimit)
+    return;
+
+  // Append new local variables corresponding to the dimensions to be converted.
+  unsigned convertCount = dimLimit - dimStart;
+  unsigned newLocalIdStart = getNumIds();
+  appendLocalId(convertCount);
+
+  // Swap the new local variables with dimensions.
+  for (unsigned i = 0; i < convertCount; ++i)
+    swapId(i + dimStart, i + newLocalIdStart);
+
+  // Remove dimensions converted to local variables.
+  removeIdRange(dimStart, dimLimit);
+}
+
+void IntegerPolyhedron::addBound(BoundType type, unsigned pos, int64_t value) {
+  assert(pos < getNumCols());
+  if (type == BoundType::EQ) {
+    unsigned row = equalities.appendExtraRow();
+    equalities(row, pos) = 1;
+    equalities(row, getNumCols() - 1) = -value;
+  } else {
+    unsigned row = inequalities.appendExtraRow();
+    inequalities(row, pos) = type == BoundType::LB ? 1 : -1;
+    inequalities(row, getNumCols() - 1) =
+        type == BoundType::LB ? -value : value;
+  }
+}
+
+void IntegerPolyhedron::addBound(BoundType type, ArrayRef<int64_t> expr,
+                                 int64_t value) {
+  assert(type != BoundType::EQ && "EQ not implemented");
+  assert(expr.size() == getNumCols());
+  unsigned row = inequalities.appendExtraRow();
+  for (unsigned i = 0, e = expr.size(); i < e; ++i)
+    inequalities(row, i) = type == BoundType::LB ? expr[i] : -expr[i];
+  inequalities(inequalities.getNumRows() - 1, getNumCols() - 1) +=
+      type == BoundType::LB ? -value : value;
+}
+
+/// Adds a new local identifier as the floordiv of an affine function of other
+/// identifiers, the coefficients of which are provided in 'dividend' and with
+/// respect to a positive constant 'divisor'. Two constraints are added to the
+/// system to capture equivalence with the floordiv.
+///      q = expr floordiv c    <=>   c*q <= expr <= c*q + c - 1.
+void IntegerPolyhedron::addLocalFloorDiv(ArrayRef<int64_t> dividend,
+                                         int64_t divisor) {
+  assert(dividend.size() == getNumCols() && "incorrect dividend size");
+  assert(divisor > 0 && "positive divisor expected");
+
+  appendLocalId();
+
+  // Add two constraints for this new identifier 'q'.
+  SmallVector<int64_t, 8> bound(dividend.size() + 1);
+
+  // dividend - q * divisor >= 0
+  std::copy(dividend.begin(), dividend.begin() + dividend.size() - 1,
+            bound.begin());
+  bound.back() = dividend.back();
+  bound[getNumIds() - 1] = -divisor;
+  addInequality(bound);
+
+  // -dividend +qdivisor * q + divisor - 1 >= 0
+  std::transform(bound.begin(), bound.end(), bound.begin(),
+                 std::negate<int64_t>());
+  bound[bound.size() - 1] += divisor - 1;
+  addInequality(bound);
+}
+
+void IntegerPolyhedron::setDimSymbolSeparation(unsigned newSymbolCount) {
+  assert(newSymbolCount <= numDims + numSymbols &&
+         "invalid separation position");
+  numDims = numDims + numSymbols - newSymbolCount;
+  numSymbols = newSymbolCount;
+}
+
+/// Finds an equality that equates the specified identifier to a constant.
+/// Returns the position of the equality row. If 'symbolic' is set to true,
+/// symbols are also treated like a constant, i.e., an affine function of the
+/// symbols is also treated like a constant. Returns -1 if such an equality
+/// could not be found.
+static int findEqualityToConstant(const IntegerPolyhedron &cst, unsigned pos,
+                                  bool symbolic = false) {
+  assert(pos < cst.getNumIds() && "invalid position");
+  for (unsigned r = 0, e = cst.getNumEqualities(); r < e; r++) {
+    int64_t v = cst.atEq(r, pos);
+    if (v * v != 1)
+      continue;
+    unsigned c;
+    unsigned f = symbolic ? cst.getNumDimIds() : cst.getNumIds();
+    // This checks for zeros in all positions other than 'pos' in [0, f)
+    for (c = 0; c < f; c++) {
+      if (c == pos)
+        continue;
+      if (cst.atEq(r, c) != 0) {
+        // Dependent on another identifier.
+        break;
+      }
+    }
+    if (c == f)
+      // Equality is free of other identifiers.
+      return r;
+  }
+  return -1;
+}
+
+void IntegerPolyhedron::setAndEliminate(unsigned pos,
+                                        ArrayRef<int64_t> values) {
+  if (values.empty())
+    return;
+  assert(pos + values.size() <= getNumIds() &&
+         "invalid position or too many values");
+  // Setting x_j = p in sum_i a_i x_i + c is equivalent to adding p*a_j to the
+  // constant term and removing the id x_j. We do this for all the ids
+  // pos, pos + 1, ... pos + values.size() - 1.
+  for (unsigned r = 0, e = getNumInequalities(); r < e; r++)
+    for (unsigned i = 0, numVals = values.size(); i < numVals; ++i)
+      atIneq(r, getNumCols() - 1) += atIneq(r, pos + i) * values[i];
+  for (unsigned r = 0, e = getNumEqualities(); r < e; r++)
+    for (unsigned i = 0, numVals = values.size(); i < numVals; ++i)
+      atEq(r, getNumCols() - 1) += atEq(r, pos + i) * values[i];
+  removeIdRange(pos, pos + values.size());
+}
+
+LogicalResult IntegerPolyhedron::constantFoldId(unsigned pos) {
+  assert(pos < getNumIds() && "invalid position");
+  int rowIdx;
+  if ((rowIdx = findEqualityToConstant(*this, pos)) == -1)
+    return failure();
+
+  // atEq(rowIdx, pos) is either -1 or 1.
+  assert(atEq(rowIdx, pos) * atEq(rowIdx, pos) == 1);
+  int64_t constVal = -atEq(rowIdx, getNumCols() - 1) / atEq(rowIdx, pos);
+  setAndEliminate(pos, constVal);
+  return success();
+}
+
+void IntegerPolyhedron::constantFoldIdRange(unsigned pos, unsigned num) {
+  for (unsigned s = pos, t = pos, e = pos + num; s < e; s++) {
+    if (failed(constantFoldId(t)))
+      t++;
+  }
+}
+
+/// Returns a non-negative constant bound on the extent (upper bound - lower
+/// bound) of the specified identifier if it is found to be a constant; returns
+/// None if it's not a constant. This methods treats symbolic identifiers
+/// specially, i.e., it looks for constant differences between affine
+/// expressions involving only the symbolic identifiers. See comments at
+/// function definition for example. 'lb', if provided, is set to the lower
+/// bound associated with the constant difference. Note that 'lb' is purely
+/// symbolic and thus will contain the coefficients of the symbolic identifiers
+/// and the constant coefficient.
+//  Egs: 0 <= i <= 15, return 16.
+//       s0 + 2 <= i <= s0 + 17, returns 16. (s0 has to be a symbol)
+//       s0 + s1 + 16 <= d0 <= s0 + s1 + 31, returns 16.
+//       s0 - 7 <= 8*j <= s0 returns 1 with lb = s0, lbDivisor = 8 (since lb =
+//       ceil(s0 - 7 / 8) = floor(s0 / 8)).
+Optional<int64_t> IntegerPolyhedron::getConstantBoundOnDimSize(
+    unsigned pos, SmallVectorImpl<int64_t> *lb, int64_t *boundFloorDivisor,
+    SmallVectorImpl<int64_t> *ub, unsigned *minLbPos,
+    unsigned *minUbPos) const {
+  assert(pos < getNumDimIds() && "Invalid identifier position");
+
+  // Find an equality for 'pos'^th identifier that equates it to some function
+  // of the symbolic identifiers (+ constant).
+  int eqPos = findEqualityToConstant(*this, pos, /*symbolic=*/true);
+  if (eqPos != -1) {
+    auto eq = getEquality(eqPos);
+    // If the equality involves a local var, punt for now.
+    // TODO: this can be handled in the future by using the explicit
+    // representation of the local vars.
+    if (!std::all_of(eq.begin() + getNumDimAndSymbolIds(), eq.end() - 1,
+                     [](int64_t coeff) { return coeff == 0; }))
+      return None;
+
+    // This identifier can only take a single value.
+    if (lb) {
+      // Set lb to that symbolic value.
+      lb->resize(getNumSymbolIds() + 1);
+      if (ub)
+        ub->resize(getNumSymbolIds() + 1);
+      for (unsigned c = 0, f = getNumSymbolIds() + 1; c < f; c++) {
+        int64_t v = atEq(eqPos, pos);
+        // atEq(eqRow, pos) is either -1 or 1.
+        assert(v * v == 1);
+        (*lb)[c] = v < 0 ? atEq(eqPos, getNumDimIds() + c) / -v
+                         : -atEq(eqPos, getNumDimIds() + c) / v;
+        // Since this is an equality, ub = lb.
+        if (ub)
+          (*ub)[c] = (*lb)[c];
+      }
+      assert(boundFloorDivisor &&
+             "both lb and divisor or none should be provided");
+      *boundFloorDivisor = 1;
+    }
+    if (minLbPos)
+      *minLbPos = eqPos;
+    if (minUbPos)
+      *minUbPos = eqPos;
+    return 1;
+  }
+
+  // Check if the identifier appears at all in any of the inequalities.
+  unsigned r, e;
+  for (r = 0, e = getNumInequalities(); r < e; r++) {
+    if (atIneq(r, pos) != 0)
+      break;
+  }
+  if (r == e)
+    // If it doesn't, there isn't a bound on it.
+    return None;
+
+  // Positions of constraints that are lower/upper bounds on the variable.
+  SmallVector<unsigned, 4> lbIndices, ubIndices;
+
+  // Gather all symbolic lower bounds and upper bounds of the variable, i.e.,
+  // the bounds can only involve symbolic (and local) identifiers. Since the
+  // canonical form c_1*x_1 + c_2*x_2 + ... + c_0 >= 0, a constraint is a lower
+  // bound for x_i if c_i >= 1, and an upper bound if c_i <= -1.
+  getLowerAndUpperBoundIndices(pos, &lbIndices, &ubIndices,
+                               /*eqIndices=*/nullptr, /*offset=*/0,
+                               /*num=*/getNumDimIds());
+
+  Optional<int64_t> minDiff = None;
+  unsigned minLbPosition = 0, minUbPosition = 0;
+  for (auto ubPos : ubIndices) {
+    for (auto lbPos : lbIndices) {
+      // Look for a lower bound and an upper bound that only differ by a
+      // constant, i.e., pairs of the form  0 <= c_pos - f(c_i's) <= diffConst.
+      // For example, if ii is the pos^th variable, we are looking for
+      // constraints like ii >= i, ii <= ii + 50, 50 being the difference. The
+      // minimum among all such constant differences is kept since that's the
+      // constant bounding the extent of the pos^th variable.
+      unsigned j, e;
+      for (j = 0, e = getNumCols() - 1; j < e; j++)
+        if (atIneq(ubPos, j) != -atIneq(lbPos, j)) {
+          break;
+        }
+      if (j < getNumCols() - 1)
+        continue;
+      int64_t diff = ceilDiv(atIneq(ubPos, getNumCols() - 1) +
+                                 atIneq(lbPos, getNumCols() - 1) + 1,
+                             atIneq(lbPos, pos));
+      // This bound is non-negative by definition.
+      diff = std::max<int64_t>(diff, 0);
+      if (minDiff == None || diff < minDiff) {
+        minDiff = diff;
+        minLbPosition = lbPos;
+        minUbPosition = ubPos;
+      }
+    }
+  }
+  if (lb && minDiff.hasValue()) {
+    // Set lb to the symbolic lower bound.
+    lb->resize(getNumSymbolIds() + 1);
+    if (ub)
+      ub->resize(getNumSymbolIds() + 1);
+    // The lower bound is the ceildiv of the lb constraint over the coefficient
+    // of the variable at 'pos'. We express the ceildiv equivalently as a floor
+    // for uniformity. For eg., if the lower bound constraint was: 32*d0 - N +
+    // 31 >= 0, the lower bound for d0 is ceil(N - 31, 32), i.e., floor(N, 32).
+    *boundFloorDivisor = atIneq(minLbPosition, pos);
+    assert(*boundFloorDivisor == -atIneq(minUbPosition, pos));
+    for (unsigned c = 0, e = getNumSymbolIds() + 1; c < e; c++) {
+      (*lb)[c] = -atIneq(minLbPosition, getNumDimIds() + c);
+    }
+    if (ub) {
+      for (unsigned c = 0, e = getNumSymbolIds() + 1; c < e; c++)
+        (*ub)[c] = atIneq(minUbPosition, getNumDimIds() + c);
+    }
+    // The lower bound leads to a ceildiv while the upper bound is a floordiv
+    // whenever the coefficient at pos != 1. ceildiv (val / d) = floordiv (val +
+    // d - 1 / d); hence, the addition of 'atIneq(minLbPosition, pos) - 1' to
+    // the constant term for the lower bound.
+    (*lb)[getNumSymbolIds()] += atIneq(minLbPosition, pos) - 1;
+  }
+  if (minLbPos)
+    *minLbPos = minLbPosition;
+  if (minUbPos)
+    *minUbPos = minUbPosition;
+  return minDiff;
+}
+
+template <bool isLower>
+Optional<int64_t>
+IntegerPolyhedron::computeConstantLowerOrUpperBound(unsigned pos) {
+  assert(pos < getNumIds() && "invalid position");
+  // Project to 'pos'.
+  projectOut(0, pos);
+  projectOut(1, getNumIds() - 1);
+  // Check if there's an equality equating the '0'^th identifier to a constant.
+  int eqRowIdx = findEqualityToConstant(*this, 0, /*symbolic=*/false);
+  if (eqRowIdx != -1)
+    // atEq(rowIdx, 0) is either -1 or 1.
+    return -atEq(eqRowIdx, getNumCols() - 1) / atEq(eqRowIdx, 0);
+
+  // Check if the identifier appears at all in any of the inequalities.
+  unsigned r, e;
+  for (r = 0, e = getNumInequalities(); r < e; r++) {
+    if (atIneq(r, 0) != 0)
+      break;
+  }
+  if (r == e)
+    // If it doesn't, there isn't a bound on it.
+    return None;
+
+  Optional<int64_t> minOrMaxConst = None;
+
+  // Take the max across all const lower bounds (or min across all constant
+  // upper bounds).
+  for (unsigned r = 0, e = getNumInequalities(); r < e; r++) {
+    if (isLower) {
+      if (atIneq(r, 0) <= 0)
+        // Not a lower bound.
+        continue;
+    } else if (atIneq(r, 0) >= 0) {
+      // Not an upper bound.
+      continue;
+    }
+    unsigned c, f;
+    for (c = 0, f = getNumCols() - 1; c < f; c++)
+      if (c != 0 && atIneq(r, c) != 0)
+        break;
+    if (c < getNumCols() - 1)
+      // Not a constant bound.
+      continue;
+
+    int64_t boundConst =
+        isLower ? mlir::ceilDiv(-atIneq(r, getNumCols() - 1), atIneq(r, 0))
+                : mlir::floorDiv(atIneq(r, getNumCols() - 1), -atIneq(r, 0));
+    if (isLower) {
+      if (minOrMaxConst == None || boundConst > minOrMaxConst)
+        minOrMaxConst = boundConst;
+    } else {
+      if (minOrMaxConst == None || boundConst < minOrMaxConst)
+        minOrMaxConst = boundConst;
+    }
+  }
+  return minOrMaxConst;
+}
+
+Optional<int64_t> IntegerPolyhedron::getConstantBound(BoundType type,
+                                                      unsigned pos) const {
+  assert(type != BoundType::EQ && "EQ not implemented");
+  IntegerPolyhedron tmpCst(*this);
+  if (type == BoundType::LB)
+    return tmpCst.computeConstantLowerOrUpperBound</*isLower=*/true>(pos);
+  return tmpCst.computeConstantLowerOrUpperBound</*isLower=*/false>(pos);
+}
+
+// A simple (naive and conservative) check for hyper-rectangularity.
+bool IntegerPolyhedron::isHyperRectangular(unsigned pos, unsigned num) const {
+  assert(pos < getNumCols() - 1);
+  // Check for two non-zero coefficients in the range [pos, pos + sum).
+  for (unsigned r = 0, e = getNumInequalities(); r < e; r++) {
+    unsigned sum = 0;
+    for (unsigned c = pos; c < pos + num; c++) {
+      if (atIneq(r, c) != 0)
+        sum++;
+    }
+    if (sum > 1)
+      return false;
+  }
+  for (unsigned r = 0, e = getNumEqualities(); r < e; r++) {
+    unsigned sum = 0;
+    for (unsigned c = pos; c < pos + num; c++) {
+      if (atEq(r, c) != 0)
+        sum++;
+    }
+    if (sum > 1)
+      return false;
+  }
+  return true;
+}
+
+/// Removes duplicate constraints, trivially true constraints, and constraints
+/// that can be detected as redundant as a result of differing only in their
+/// constant term part. A constraint of the form <non-negative constant> >= 0 is
+/// considered trivially true.
+//  Uses a DenseSet to hash and detect duplicates followed by a linear scan to
+//  remove duplicates in place.
+void IntegerPolyhedron::removeTrivialRedundancy() {
+  gcdTightenInequalities();
+  normalizeConstraintsByGCD();
+
+  // A map used to detect redundancy stemming from constraints that only differ
+  // in their constant term. The value stored is <row position, const term>
+  // for a given row.
+  SmallDenseMap<ArrayRef<int64_t>, std::pair<unsigned, int64_t>>
+      rowsWithoutConstTerm;
+  // To unique rows.
+  SmallDenseSet<ArrayRef<int64_t>, 8> rowSet;
+
+  // Check if constraint is of the form <non-negative-constant> >= 0.
+  auto isTriviallyValid = [&](unsigned r) -> bool {
+    for (unsigned c = 0, e = getNumCols() - 1; c < e; c++) {
+      if (atIneq(r, c) != 0)
+        return false;
+    }
+    return atIneq(r, getNumCols() - 1) >= 0;
+  };
+
+  // Detect and mark redundant constraints.
+  SmallVector<bool, 256> redunIneq(getNumInequalities(), false);
+  for (unsigned r = 0, e = getNumInequalities(); r < e; r++) {
+    int64_t *rowStart = &inequalities(r, 0);
+    auto row = ArrayRef<int64_t>(rowStart, getNumCols());
+    if (isTriviallyValid(r) || !rowSet.insert(row).second) {
+      redunIneq[r] = true;
+      continue;
+    }
+
+    // Among constraints that only differ in the constant term part, mark
+    // everything other than the one with the smallest constant term redundant.
+    // (eg: among i - 16j - 5 >= 0, i - 16j - 1 >=0, i - 16j - 7 >= 0, the
+    // former two are redundant).
+    int64_t constTerm = atIneq(r, getNumCols() - 1);
+    auto rowWithoutConstTerm = ArrayRef<int64_t>(rowStart, getNumCols() - 1);
+    const auto &ret =
+        rowsWithoutConstTerm.insert({rowWithoutConstTerm, {r, constTerm}});
+    if (!ret.second) {
+      // Check if the other constraint has a higher constant term.
+      auto &val = ret.first->second;
+      if (val.second > constTerm) {
+        // The stored row is redundant. Mark it so, and update with this one.
+        redunIneq[val.first] = true;
+        val = {r, constTerm};
+      } else {
+        // The one stored makes this one redundant.
+        redunIneq[r] = true;
+      }
+    }
+  }
+
+  // Scan to get rid of all rows marked redundant, in-place.
+  unsigned pos = 0;
+  for (unsigned r = 0, e = getNumInequalities(); r < e; r++)
+    if (!redunIneq[r])
+      inequalities.copyRow(r, pos++);
+
+  inequalities.resizeVertically(pos);
+
+  // TODO: consider doing this for equalities as well, but probably not worth
+  // the savings.
+}
+
+static std::pair<unsigned, unsigned>
+getNewNumDimsSymbols(unsigned pos, const IntegerPolyhedron &cst) {
+  unsigned numDims = cst.getNumDimIds();
+  unsigned numSymbols = cst.getNumSymbolIds();
+  unsigned newNumDims, newNumSymbols;
+  if (pos < numDims) {
+    newNumDims = numDims - 1;
+    newNumSymbols = numSymbols;
+  } else if (pos < numDims + numSymbols) {
+    assert(numSymbols >= 1);
+    newNumDims = numDims;
+    newNumSymbols = numSymbols - 1;
+  } else {
+    newNumDims = numDims;
+    newNumSymbols = numSymbols;
+  }
+  return {newNumDims, newNumSymbols};
+}
+
+void IntegerPolyhedron::clearAndCopyFrom(const IntegerPolyhedron &other) {
+  *this = other;
+}
+
+#undef DEBUG_TYPE
+#define DEBUG_TYPE "fm"
+
+/// Eliminates identifier at the specified position using Fourier-Motzkin
+/// variable elimination. This technique is exact for rational spaces but
+/// conservative (in "rare" cases) for integer spaces. The operation corresponds
+/// to a projection operation yielding the (convex) set of integer points
+/// contained in the rational shadow of the set. An emptiness test that relies
+/// on this method will guarantee emptiness, i.e., it disproves the existence of
+/// a solution if it says it's empty.
+/// If a non-null isResultIntegerExact is passed, it is set to true if the
+/// result is also integer exact. If it's set to false, the obtained solution
+/// *may* not be exact, i.e., it may contain integer points that do not have an
+/// integer pre-image in the original set.
+///
+/// Eg:
+/// j >= 0, j <= i + 1
+/// i >= 0, i <= N + 1
+/// Eliminating i yields,
+///   j >= 0, 0 <= N + 1, j - 1 <= N + 1
+///
+/// If darkShadow = true, this method computes the dark shadow on elimination;
+/// the dark shadow is a convex integer subset of the exact integer shadow. A
+/// non-empty dark shadow proves the existence of an integer solution. The
+/// elimination in such a case could however be an under-approximation, and thus
+/// should not be used for scanning sets or used by itself for dependence
+/// checking.
+///
+/// Eg: 2-d set, * represents grid points, 'o' represents a point in the set.
+///            ^
+///            |
+///            | * * * * o o
+///         i  | * * o o o o
+///            | o * * * * *
+///            --------------->
+///                 j ->
+///
+/// Eliminating i from this system (projecting on the j dimension):
+/// rational shadow / integer light shadow:  1 <= j <= 6
+/// dark shadow:                             3 <= j <= 6
+/// exact integer shadow:                    j = 1 \union  3 <= j <= 6
+/// holes/splinters:                         j = 2
+///
+/// darkShadow = false, isResultIntegerExact = nullptr are default values.
+// TODO: a slight modification to yield dark shadow version of FM (tightened),
+// which can prove the existence of a solution if there is one.
+void IntegerPolyhedron::fourierMotzkinEliminate(unsigned pos, bool darkShadow,
+                                                bool *isResultIntegerExact) {
+  LLVM_DEBUG(llvm::dbgs() << "FM input (eliminate pos " << pos << "):\n");
+  LLVM_DEBUG(dump());
+  assert(pos < getNumIds() && "invalid position");
+  assert(hasConsistentState());
+
+  // Check if this identifier can be eliminated through a substitution.
+  for (unsigned r = 0, e = getNumEqualities(); r < e; r++) {
+    if (atEq(r, pos) != 0) {
+      // Use Gaussian elimination here (since we have an equality).
+      LogicalResult ret = gaussianEliminateId(pos);
+      (void)ret;
+      assert(succeeded(ret) && "Gaussian elimination guaranteed to succeed");
+      LLVM_DEBUG(llvm::dbgs() << "FM output (through Gaussian elimination):\n");
+      LLVM_DEBUG(dump());
+      return;
+    }
+  }
+
+  // A fast linear time tightening.
+  gcdTightenInequalities();
+
+  // Check if the identifier appears at all in any of the inequalities.
+  unsigned r, e;
+  for (r = 0, e = getNumInequalities(); r < e; r++) {
+    if (atIneq(r, pos) != 0)
+      break;
+  }
+  if (r == getNumInequalities()) {
+    // If it doesn't appear, just remove the column and return.
+    // TODO: refactor removeColumns to use it from here.
+    removeId(pos);
+    LLVM_DEBUG(llvm::dbgs() << "FM output:\n");
+    LLVM_DEBUG(dump());
+    return;
+  }
+
+  // Positions of constraints that are lower bounds on the variable.
+  SmallVector<unsigned, 4> lbIndices;
+  // Positions of constraints that are lower bounds on the variable.
+  SmallVector<unsigned, 4> ubIndices;
+  // Positions of constraints that do not involve the variable.
+  std::vector<unsigned> nbIndices;
+  nbIndices.reserve(getNumInequalities());
+
+  // Gather all lower bounds and upper bounds of the variable. Since the
+  // canonical form c_1*x_1 + c_2*x_2 + ... + c_0 >= 0, a constraint is a lower
+  // bound for x_i if c_i >= 1, and an upper bound if c_i <= -1.
+  for (unsigned r = 0, e = getNumInequalities(); r < e; r++) {
+    if (atIneq(r, pos) == 0) {
+      // Id does not appear in bound.
+      nbIndices.push_back(r);
+    } else if (atIneq(r, pos) >= 1) {
+      // Lower bound.
+      lbIndices.push_back(r);
+    } else {
+      // Upper bound.
+      ubIndices.push_back(r);
+    }
+  }
+
+  // Set the number of dimensions, symbols in the resulting system.
+  const auto &dimsSymbols = getNewNumDimsSymbols(pos, *this);
+  unsigned newNumDims = dimsSymbols.first;
+  unsigned newNumSymbols = dimsSymbols.second;
+
+  /// Create the new system which has one identifier less.
+  IntegerPolyhedron newFac(
+      lbIndices.size() * ubIndices.size() + nbIndices.size(),
+      getNumEqualities(), getNumCols() - 1, newNumDims, newNumSymbols,
+      /*numLocals=*/getNumIds() - 1 - newNumDims - newNumSymbols);
+
+  // This will be used to check if the elimination was integer exact.
+  unsigned lcmProducts = 1;
+
+  // Let x be the variable we are eliminating.
+  // For each lower bound, lb <= c_l*x, and each upper bound c_u*x <= ub, (note
+  // that c_l, c_u >= 1) we have:
+  // lb*lcm(c_l, c_u)/c_l <= lcm(c_l, c_u)*x <= ub*lcm(c_l, c_u)/c_u
+  // We thus generate a constraint:
+  // lcm(c_l, c_u)/c_l*lb <= lcm(c_l, c_u)/c_u*ub.
+  // Note if c_l = c_u = 1, all integer points captured by the resulting
+  // constraint correspond to integer points in the original system (i.e., they
+  // have integer pre-images). Hence, if the lcm's are all 1, the elimination is
+  // integer exact.
+  for (auto ubPos : ubIndices) {
+    for (auto lbPos : lbIndices) {
+      SmallVector<int64_t, 4> ineq;
+      ineq.reserve(newFac.getNumCols());
+      int64_t lbCoeff = atIneq(lbPos, pos);
+      // Note that in the comments above, ubCoeff is the negation of the
+      // coefficient in the canonical form as the view taken here is that of the
+      // term being moved to the other size of '>='.
+      int64_t ubCoeff = -atIneq(ubPos, pos);
+      // TODO: refactor this loop to avoid all branches inside.
+      for (unsigned l = 0, e = getNumCols(); l < e; l++) {
+        if (l == pos)
+          continue;
+        assert(lbCoeff >= 1 && ubCoeff >= 1 && "bounds wrongly identified");
+        int64_t lcm = mlir::lcm(lbCoeff, ubCoeff);
+        ineq.push_back(atIneq(ubPos, l) * (lcm / ubCoeff) +
+                       atIneq(lbPos, l) * (lcm / lbCoeff));
+        lcmProducts *= lcm;
+      }
+      if (darkShadow) {
+        // The dark shadow is a convex subset of the exact integer shadow. If
+        // there is a point here, it proves the existence of a solution.
+        ineq[ineq.size() - 1] += lbCoeff * ubCoeff - lbCoeff - ubCoeff + 1;
+      }
+      // TODO: we need to have a way to add inequalities in-place in
+      // IntegerPolyhedron instead of creating and copying over.
+      newFac.addInequality(ineq);
+    }
+  }
+
+  LLVM_DEBUG(llvm::dbgs() << "FM isResultIntegerExact: " << (lcmProducts == 1)
+                          << "\n");
+  if (lcmProducts == 1 && isResultIntegerExact)
+    *isResultIntegerExact = true;
+
+  // Copy over the constraints not involving this variable.
+  for (auto nbPos : nbIndices) {
+    SmallVector<int64_t, 4> ineq;
+    ineq.reserve(getNumCols() - 1);
+    for (unsigned l = 0, e = getNumCols(); l < e; l++) {
+      if (l == pos)
+        continue;
+      ineq.push_back(atIneq(nbPos, l));
+    }
+    newFac.addInequality(ineq);
+  }
+
+  assert(newFac.getNumConstraints() ==
+         lbIndices.size() * ubIndices.size() + nbIndices.size());
+
+  // Copy over the equalities.
+  for (unsigned r = 0, e = getNumEqualities(); r < e; r++) {
+    SmallVector<int64_t, 4> eq;
+    eq.reserve(newFac.getNumCols());
+    for (unsigned l = 0, e = getNumCols(); l < e; l++) {
+      if (l == pos)
+        continue;
+      eq.push_back(atEq(r, l));
+    }
+    newFac.addEquality(eq);
+  }
+
+  // GCD tightening and normalization allows detection of more trivially
+  // redundant constraints.
+  newFac.gcdTightenInequalities();
+  newFac.normalizeConstraintsByGCD();
+  newFac.removeTrivialRedundancy();
+  clearAndCopyFrom(newFac);
+  LLVM_DEBUG(llvm::dbgs() << "FM output:\n");
+  LLVM_DEBUG(dump());
+}
+
+#undef DEBUG_TYPE
+#define DEBUG_TYPE "presburger"
+
+void IntegerPolyhedron::projectOut(unsigned pos, unsigned num) {
+  if (num == 0)
+    return;
+
+  // 'pos' can be at most getNumCols() - 2 if num > 0.
+  assert((getNumCols() < 2 || pos <= getNumCols() - 2) && "invalid position");
+  assert(pos + num < getNumCols() && "invalid range");
+
+  // Eliminate as many identifiers as possible using Gaussian elimination.
+  unsigned currentPos = pos;
+  unsigned numToEliminate = num;
+  unsigned numGaussianEliminated = 0;
+
+  while (currentPos < getNumIds()) {
+    unsigned curNumEliminated =
+        gaussianEliminateIds(currentPos, currentPos + numToEliminate);
+    ++currentPos;
+    numToEliminate -= curNumEliminated + 1;
+    numGaussianEliminated += curNumEliminated;
+  }
+
+  // Eliminate the remaining using Fourier-Motzkin.
+  for (unsigned i = 0; i < num - numGaussianEliminated; i++) {
+    unsigned numToEliminate = num - numGaussianEliminated - i;
+    fourierMotzkinEliminate(
+        getBestIdToEliminate(*this, pos, pos + numToEliminate));
+  }
+
+  // Fast/trivial simplifications.
+  gcdTightenInequalities();
+  // Normalize constraints after tightening since the latter impacts this, but
+  // not the other way round.
+  normalizeConstraintsByGCD();
+}
+
+namespace {
+
+enum BoundCmpResult { Greater, Less, Equal, Unknown };
+
+/// Compares two affine bounds whose coefficients are provided in 'first' and
+/// 'second'. The last coefficient is the constant term.
+static BoundCmpResult compareBounds(ArrayRef<int64_t> a, ArrayRef<int64_t> b) {
+  assert(a.size() == b.size());
+
+  // For the bounds to be comparable, their corresponding identifier
+  // coefficients should be equal; the constant terms are then compared to
+  // determine less/greater/equal.
+
+  if (!std::equal(a.begin(), a.end() - 1, b.begin()))
+    return Unknown;
+
+  if (a.back() == b.back())
+    return Equal;
+
+  return a.back() < b.back() ? Less : Greater;
+}
+} // namespace
+
+// Returns constraints that are common to both A & B.
+static void getCommonConstraints(const IntegerPolyhedron &a,
+                                 const IntegerPolyhedron &b,
+                                 IntegerPolyhedron &c) {
+  c.reset(a.getNumDimIds(), a.getNumSymbolIds(), a.getNumLocalIds());
+  // a naive O(n^2) check should be enough here given the input sizes.
+  for (unsigned r = 0, e = a.getNumInequalities(); r < e; ++r) {
+    for (unsigned s = 0, f = b.getNumInequalities(); s < f; ++s) {
+      if (a.getInequality(r) == b.getInequality(s)) {
+        c.addInequality(a.getInequality(r));
+        break;
+      }
+    }
+  }
+  for (unsigned r = 0, e = a.getNumEqualities(); r < e; ++r) {
+    for (unsigned s = 0, f = b.getNumEqualities(); s < f; ++s) {
+      if (a.getEquality(r) == b.getEquality(s)) {
+        c.addEquality(a.getEquality(r));
+        break;
+      }
+    }
+  }
+}
+
+// Computes the bounding box with respect to 'other' by finding the min of the
+// lower bounds and the max of the upper bounds along each of the dimensions.
+LogicalResult
+IntegerPolyhedron::unionBoundingBox(const IntegerPolyhedron &otherCst) {
+  assert(otherCst.getNumDimIds() == numDims && "dims mismatch");
+  assert(otherCst.getNumLocalIds() == 0 && "local ids not supported here");
+  assert(getNumLocalIds() == 0 && "local ids not supported yet here");
+
+  // Get the constraints common to both systems; these will be added as is to
+  // the union.
+  IntegerPolyhedron commonCst;
+  getCommonConstraints(*this, otherCst, commonCst);
+
+  std::vector<SmallVector<int64_t, 8>> boundingLbs;
+  std::vector<SmallVector<int64_t, 8>> boundingUbs;
+  boundingLbs.reserve(2 * getNumDimIds());
+  boundingUbs.reserve(2 * getNumDimIds());
+
+  // To hold lower and upper bounds for each dimension.
+  SmallVector<int64_t, 4> lb, otherLb, ub, otherUb;
+  // To compute min of lower bounds and max of upper bounds for each dimension.
+  SmallVector<int64_t, 4> minLb(getNumSymbolIds() + 1);
+  SmallVector<int64_t, 4> maxUb(getNumSymbolIds() + 1);
+  // To compute final new lower and upper bounds for the union.
+  SmallVector<int64_t, 8> newLb(getNumCols()), newUb(getNumCols());
+
+  int64_t lbFloorDivisor, otherLbFloorDivisor;
+  for (unsigned d = 0, e = getNumDimIds(); d < e; ++d) {
+    auto extent = getConstantBoundOnDimSize(d, &lb, &lbFloorDivisor, &ub);
+    if (!extent.hasValue())
+      // TODO: symbolic extents when necessary.
+      // TODO: handle union if a dimension is unbounded.
+      return failure();
+
+    auto otherExtent = otherCst.getConstantBoundOnDimSize(
+        d, &otherLb, &otherLbFloorDivisor, &otherUb);
+    if (!otherExtent.hasValue() || lbFloorDivisor != otherLbFloorDivisor)
+      // TODO: symbolic extents when necessary.
+      return failure();
+
+    assert(lbFloorDivisor > 0 && "divisor always expected to be positive");
+
+    auto res = compareBounds(lb, otherLb);
+    // Identify min.
+    if (res == BoundCmpResult::Less || res == BoundCmpResult::Equal) {
+      minLb = lb;
+      // Since the divisor is for a floordiv, we need to convert to ceildiv,
+      // i.e., i >= expr floordiv div <=> i >= (expr - div + 1) ceildiv div <=>
+      // div * i >= expr - div + 1.
+      minLb.back() -= lbFloorDivisor - 1;
+    } else if (res == BoundCmpResult::Greater) {
+      minLb = otherLb;
+      minLb.back() -= otherLbFloorDivisor - 1;
+    } else {
+      // Uncomparable - check for constant lower/upper bounds.
+      auto constLb = getConstantBound(BoundType::LB, d);
+      auto constOtherLb = otherCst.getConstantBound(BoundType::LB, d);
+      if (!constLb.hasValue() || !constOtherLb.hasValue())
+        return failure();
+      std::fill(minLb.begin(), minLb.end(), 0);
+      minLb.back() = std::min(constLb.getValue(), constOtherLb.getValue());
+    }
+
+    // Do the same for ub's but max of upper bounds. Identify max.
+    auto uRes = compareBounds(ub, otherUb);
+    if (uRes == BoundCmpResult::Greater || uRes == BoundCmpResult::Equal) {
+      maxUb = ub;
+    } else if (uRes == BoundCmpResult::Less) {
+      maxUb = otherUb;
+    } else {
+      // Uncomparable - check for constant lower/upper bounds.
+      auto constUb = getConstantBound(BoundType::UB, d);
+      auto constOtherUb = otherCst.getConstantBound(BoundType::UB, d);
+      if (!constUb.hasValue() || !constOtherUb.hasValue())
+        return failure();
+      std::fill(maxUb.begin(), maxUb.end(), 0);
+      maxUb.back() = std::max(constUb.getValue(), constOtherUb.getValue());
+    }
+
+    std::fill(newLb.begin(), newLb.end(), 0);
+    std::fill(newUb.begin(), newUb.end(), 0);
+
+    // The divisor for lb, ub, otherLb, otherUb at this point is lbDivisor,
+    // and so it's the divisor for newLb and newUb as well.
+    newLb[d] = lbFloorDivisor;
+    newUb[d] = -lbFloorDivisor;
+    // Copy over the symbolic part + constant term.
+    std::copy(minLb.begin(), minLb.end(), newLb.begin() + getNumDimIds());
+    std::transform(newLb.begin() + getNumDimIds(), newLb.end(),
+                   newLb.begin() + getNumDimIds(), std::negate<int64_t>());
+    std::copy(maxUb.begin(), maxUb.end(), newUb.begin() + getNumDimIds());
+
+    boundingLbs.push_back(newLb);
+    boundingUbs.push_back(newUb);
+  }
+
+  // Clear all constraints and add the lower/upper bounds for the bounding box.
+  clearConstraints();
+  for (unsigned d = 0, e = getNumDimIds(); d < e; ++d) {
+    addInequality(boundingLbs[d]);
+    addInequality(boundingUbs[d]);
+  }
+
+  // Add the constraints that were common to both systems.
+  append(commonCst);
+  removeTrivialRedundancy();
+
+  // TODO: copy over pure symbolic constraints from this and 'other' over to the
+  // union (since the above are just the union along dimensions); we shouldn't
+  // be discarding any other constraints on the symbols.
+
+  return success();
+}
+
+bool IntegerPolyhedron::isColZero(unsigned pos) const {
+  unsigned rowPos;
+  return !findConstraintWithNonZeroAt(pos, /*isEq=*/false, &rowPos) &&
+         !findConstraintWithNonZeroAt(pos, /*isEq=*/true, &rowPos);
+}
+
+/// Find positions of inequalities and equalities that do not have a coefficient
+/// for [pos, pos + num) identifiers.
+static void getIndependentConstraints(const IntegerPolyhedron &cst,
+                                      unsigned pos, unsigned num,
+                                      SmallVectorImpl<unsigned> &nbIneqIndices,
+                                      SmallVectorImpl<unsigned> &nbEqIndices) {
+  assert(pos < cst.getNumIds() && "invalid start position");
+  assert(pos + num <= cst.getNumIds() && "invalid limit");
+
+  for (unsigned r = 0, e = cst.getNumInequalities(); r < e; r++) {
+    // The bounds are to be independent of [offset, offset + num) columns.
+    unsigned c;
+    for (c = pos; c < pos + num; ++c) {
+      if (cst.atIneq(r, c) != 0)
+        break;
+    }
+    if (c == pos + num)
+      nbIneqIndices.push_back(r);
+  }
+
+  for (unsigned r = 0, e = cst.getNumEqualities(); r < e; r++) {
+    // The bounds are to be independent of [offset, offset + num) columns.
+    unsigned c;
+    for (c = pos; c < pos + num; ++c) {
+      if (cst.atEq(r, c) != 0)
+        break;
+    }
+    if (c == pos + num)
+      nbEqIndices.push_back(r);
+  }
+}
+
+void IntegerPolyhedron::removeIndependentConstraints(unsigned pos,
+                                                     unsigned num) {
+  assert(pos + num <= getNumIds() && "invalid range");
+
+  // Remove constraints that are independent of these identifiers.
+  SmallVector<unsigned, 4> nbIneqIndices, nbEqIndices;
+  getIndependentConstraints(*this, /*pos=*/0, num, nbIneqIndices, nbEqIndices);
+
+  // Iterate in reverse so that indices don't have to be updated.
+  // TODO: This method can be made more efficient (because removal of each
+  // inequality leads to much shifting/copying in the underlying buffer).
+  for (auto nbIndex : llvm::reverse(nbIneqIndices))
+    removeInequality(nbIndex);
+  for (auto nbIndex : llvm::reverse(nbEqIndices))
+    removeEquality(nbIndex);
+}
+
+void IntegerPolyhedron::printSpace(raw_ostream &os) const {
+  os << "\nConstraints (" << getNumDimIds() << " dims, " << getNumSymbolIds()
+     << " symbols, " << getNumLocalIds() << " locals), (" << getNumConstraints()
+     << " constraints)\n";
+  os << "(";
+  os << " const)\n";
+}
+
+void IntegerPolyhedron::print(raw_ostream &os) const {
+  assert(hasConsistentState());
+  printSpace(os);
+  for (unsigned i = 0, e = getNumEqualities(); i < e; ++i) {
+    for (unsigned j = 0, f = getNumCols(); j < f; ++j) {
+      os << atEq(i, j) << " ";
+    }
+    os << "= 0\n";
+  }
+  for (unsigned i = 0, e = getNumInequalities(); i < e; ++i) {
+    for (unsigned j = 0, f = getNumCols(); j < f; ++j) {
+      os << atIneq(i, j) << " ";
+    }
+    os << ">= 0\n";
+  }
+  os << '\n';
+}
+
+void IntegerPolyhedron::dump() const { print(llvm::errs()); }

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/LinearTransform.h"
-#include "mlir/Analysis/AffineStructures.h"
+#include "mlir/Analysis/Presburger/LinearTransform.h"
+#include "mlir/Analysis/Presburger/IntegerPolyhedron.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
LinearTransform::preMultiplyColumn(ArrayRef<int64_t> colVec) const {
 
   SmallVector<int64_t, 8> result(matrix.getNumRows(), 0);
   for (unsigned row = 0, e = matrix.getNumRows(); row < e; row++)
     for (unsigned i = 0, e = matrix.getNumColumns(); i < e; i++)
       result[row] += matrix(row, i) * colVec[i];
   return result;
 }
 
-FlatAffineConstraints
-LinearTransform::applyTo(const FlatAffineConstraints &fac) const {
-  FlatAffineConstraints result(fac.getNumIds());
+IntegerPolyhedron LinearTransform::applyTo(const IntegerPolyhedron &fac) const {
+  IntegerPolyhedron result(fac.getNumIds());
 
   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 = postMultiplyRow(eq.drop_back());
     newEq.push_back(c);

mlir/lib/Analysis/Presburger/Simplex.cpp

SimplexBase::SimplexBase(unsigned nVar)
   colUnknown.push_back(nullIndex);
   for (unsigned i = 0; i < nVar; ++i) {
     var.emplace_back(Orientation::Column, /*restricted=*/false, /*pos=*/nCol);
     colUnknown.push_back(i);
     nCol++;
   }
 }
 
-SimplexBase::SimplexBase(const FlatAffineConstraints &constraints)
+SimplexBase::SimplexBase(const IntegerPolyhedron &constraints)
     : SimplexBase(constraints.getNumIds()) {
   for (unsigned i = 0, numIneqs = constraints.getNumInequalities();
        i < numIneqs; ++i)
     addInequality(constraints.getInequality(i));
   for (unsigned i = 0, numEqs = constraints.getNumEqualities(); i < numEqs; ++i)
     addEquality(constraints.getEquality(i));
 }
 
for (unsigned i = 0; i < count; ++i) {
     var.emplace_back(Orientation::Column, /*restricted=*/false,
                      /*pos=*/nCol - 1);
     colUnknown.push_back(var.size() - 1);
   }
   tableau.resizeHorizontally(nCol);
   undoLog.insert(undoLog.end(), count, UndoLogEntry::RemoveLastVariable);
 }
 
-/// Add all the constraints from the given FlatAffineConstraints.
-void SimplexBase::intersectFlatAffineConstraints(
-    const FlatAffineConstraints &fac) {
+/// Add all the constraints from the given IntegerPolyhedron.
+void SimplexBase::intersectIntegerPolyhedron(const IntegerPolyhedron &fac) {
   assert(fac.getNumIds() == getNumVariables() &&
-         "FlatAffineConstraints must have same dimensionality as simplex");
+         "IntegerPolyhedron must have same dimensionality as simplex");
   for (unsigned i = 0, e = fac.getNumInequalities(); i < e; ++i)
     addInequality(fac.getInequality(i));
   for (unsigned i = 0, e = fac.getNumEqualities(); i < e; ++i)
     addEquality(fac.getEquality(i));
 }
 
 Optional<Fraction> Simplex::computeRowOptimum(Direction direction,
                                               unsigned row) {
for (unsigned col = 0; col < nCol; ++col)
       os << tableau(row, col) << '\t';
     os << '\n';
   }
   os << '\n';
 }
 
 void SimplexBase::dump() const { print(llvm::errs()); }
 
-bool Simplex::isRationalSubsetOf(const FlatAffineConstraints &fac) {
+bool Simplex::isRationalSubsetOf(const IntegerPolyhedron &fac) {
   if (isEmpty())
     return true;
 
   for (unsigned i = 0, e = fac.getNumInequalities(); i < e; ++i)
     if (!isRedundantInequality(fac.getInequality(i)))
       return false;
 
   for (unsigned i = 0, e = fac.getNumEqualities(); i < e; ++i)

mlir/lib/Analysis/PresburgerSet.cpp

 //===- Set.cpp - MLIR PresburgerSet 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/PresburgerSet.h"
 #include "mlir/Analysis/Presburger/Simplex.h"
 #include "llvm/ADT/STLExtras.h"
 #include "llvm/ADT/SmallBitVector.h"
 
 using namespace mlir;
 
-PresburgerSet::PresburgerSet(const FlatAffineConstraints &fac)
+PresburgerSet::PresburgerSet(const IntegerPolyhedron &fac)
     : nDim(fac.getNumDimIds()), nSym(fac.getNumSymbolIds()) {
   unionFACInPlace(fac);
 }
 
 unsigned PresburgerSet::getNumFACs() const {
   return flatAffineConstraints.size();
 }
 
 unsigned PresburgerSet::getNumDims() const { return nDim; }
 
 unsigned PresburgerSet::getNumSyms() const { return nSym; }
 
-ArrayRef<FlatAffineConstraints>
-PresburgerSet::getAllFlatAffineConstraints() const {
+ArrayRef<IntegerPolyhedron> PresburgerSet::getAllIntegerPolyhedron() const {
   return flatAffineConstraints;
 }
 
-const FlatAffineConstraints &
-PresburgerSet::getFlatAffineConstraints(unsigned index) const {
+const IntegerPolyhedron &
+PresburgerSet::getIntegerPolyhedron(unsigned index) const {
   assert(index < flatAffineConstraints.size() && "index out of bounds!");
   return flatAffineConstraints[index];
 }
 
-/// Assert that the FlatAffineConstraints and PresburgerSet live in
+/// Assert that the IntegerPolyhedron and PresburgerSet live in
 /// compatible spaces.
-static void assertDimensionsCompatible(const FlatAffineConstraints &fac,
+static void assertDimensionsCompatible(const IntegerPolyhedron &fac,
                                        const PresburgerSet &set) {
   assert(fac.getNumDimIds() == set.getNumDims() &&
-         "Number of dimensions of the FlatAffineConstraints and PresburgerSet"
+         "Number of dimensions of the IntegerPolyhedron and PresburgerSet"
          "do not match!");
   assert(fac.getNumSymbolIds() == set.getNumSyms() &&
-         "Number of symbols of the FlatAffineConstraints and PresburgerSet"
+         "Number of symbols of the IntegerPolyhedron and PresburgerSet"
          "do not match!");
 }
 
 /// Assert that the two PresburgerSets live in compatible spaces.
 static void assertDimensionsCompatible(const PresburgerSet &setA,
                                        const PresburgerSet &setB) {
   assert(setA.getNumDims() == setB.getNumDims() &&
          "Number of dimensions of the PresburgerSets do not match!");
   assert(setA.getNumSyms() == setB.getNumSyms() &&
          "Number of symbols of the PresburgerSets do not match!");
 }
 
 /// Mutate this set, turning it into the union of this set and the given
-/// FlatAffineConstraints.
-void PresburgerSet::unionFACInPlace(const FlatAffineConstraints &fac) {
+/// IntegerPolyhedron.
+void PresburgerSet::unionFACInPlace(const IntegerPolyhedron &fac) {
   assertDimensionsCompatible(fac, *this);
   flatAffineConstraints.push_back(fac);
 }
 
 /// Mutate this set, turning it into the union of this set and the given set.
 ///
 /// This is accomplished by simply adding all the FACs of the given set to this
 /// set.
 void PresburgerSet::unionSetInPlace(const PresburgerSet &set) {
   assertDimensionsCompatible(set, *this);
-  for (const FlatAffineConstraints &fac : set.flatAffineConstraints)
+  for (const IntegerPolyhedron &fac : set.flatAffineConstraints)
     unionFACInPlace(fac);
 }
 
 /// Return the union of this set and the given set.
 PresburgerSet PresburgerSet::unionSet(const PresburgerSet &set) const {
   assertDimensionsCompatible(set, *this);
   PresburgerSet result = *this;
   result.unionSetInPlace(set);
   return result;
 }
 
 /// A point is contained in the union iff any of the parts contain the point.
 bool PresburgerSet::containsPoint(ArrayRef<int64_t> point) const {
-  for (const FlatAffineConstraints &fac : flatAffineConstraints) {
+  for (const IntegerPolyhedron &fac : flatAffineConstraints) {
     if (fac.containsPoint(point))
       return true;
   }
   return false;
 }
 
 PresburgerSet PresburgerSet::getUniverse(unsigned nDim, unsigned nSym) {
   PresburgerSet result(nDim, nSym);
-  result.unionFACInPlace(FlatAffineConstraints::getUniverse(nDim, nSym));
+  result.unionFACInPlace(IntegerPolyhedron::getUniverse(nDim, nSym));
   return result;
 }
 
 PresburgerSet PresburgerSet::getEmptySet(unsigned nDim, unsigned nSym) {
   return PresburgerSet(nDim, nSym);
 }
 
 // Return the intersection of this set with the given set.
 //
 // We directly compute (S_1 or S_2 ...) and (T_1 or T_2 ...)
 // as (S_1 and T_1) or (S_1 and T_2) or ...
 //
 // If S_i or T_j have local variables, then S_i and T_j contains the local
 // variables of both.
 PresburgerSet PresburgerSet::intersect(const PresburgerSet &set) const {
   assertDimensionsCompatible(set, *this);
 
   PresburgerSet result(nDim, nSym);
-  for (const FlatAffineConstraints &csA : flatAffineConstraints) {
-    for (const FlatAffineConstraints &csB : set.flatAffineConstraints) {
-      FlatAffineConstraints csACopy = csA, csBCopy = csB;
+  for (const IntegerPolyhedron &csA : flatAffineConstraints) {
+    for (const IntegerPolyhedron &csB : set.flatAffineConstraints) {
+      IntegerPolyhedron csACopy = csA, csBCopy = csB;
       csACopy.mergeLocalIds(csBCopy);
       csACopy.append(std::move(csBCopy));
       if (!csACopy.isEmpty())
         result.unionFACInPlace(std::move(csACopy));
     }
   }
   return result;
 }
static SmallVector<int64_t, 8> getComplementIneq(ArrayRef<int64_t> ineq) {
   return coeffs;
 }
 
 /// Return the set difference b \ s and accumulate the result into `result`.
 /// `simplex` must correspond to b.
 ///
 /// In the following, U denotes union, ^ denotes intersection, \ denotes set
 /// difference and ~ denotes complement.
-/// Let b be the FlatAffineConstraints and s = (U_i s_i) be the set. We want
+/// Let b be the IntegerPolyhedron and s = (U_i s_i) be the set. We want
 /// b \ (U_i s_i).
 ///
 /// Let s_i = ^_j s_ij, where each s_ij is a single inequality. To compute
 /// b \ s_i = b ^ ~s_i, we partition s_i based on the first violated inequality:
 /// ~s_i = (~s_i1) U (s_i1 ^ ~s_i2) U (s_i1 ^ s_i2 ^ ~s_i3) U ...
 /// And the required result is (b ^ ~s_i1) U (b ^ s_i1 ^ ~s_i2) U ...
 /// We recurse by subtracting U_{j > i} S_j from each of these parts and
 /// returning the union of the results. Each equality is handled as a

 ///
 /// As a heuristic, we try adding all the constraints and check if simplex
 /// says that the intersection is empty. If it is, then subtracting this FAC is
 /// a no-op and we just skip it. Also, in the process we find out that some
 /// constraints are redundant. These redundant constraints are ignored.
 ///
 /// b and simplex are callee saved, i.e., their values on return are
 /// semantically equivalent to their values when the function is called.
-static void subtractRecursively(FlatAffineConstraints &b, Simplex &simplex,
+static void subtractRecursively(IntegerPolyhedron &b, Simplex &simplex,
                                 const PresburgerSet &s, unsigned i,
                                 PresburgerSet &result) {
   if (i == s.getNumFACs()) {
     result.unionFACInPlace(b);
     return;
   }
-  FlatAffineConstraints sI = s.getFlatAffineConstraints(i);
+  IntegerPolyhedron sI = s.getIntegerPolyhedron(i);
 
   // Below, we append some additional constraints and ids to b. We want to
   // rollback b to its initial state before returning, which we will do by
   // removing all constraints beyond the original number of inequalities
   // and equalities, so we store these counts first.
   const unsigned bInitNumIneqs = b.getNumInequalities();
   const unsigned bInitNumEqs = b.getNumEqualities();
   const unsigned bInitNumLocals = b.getNumLocalIds();
   // Similarly, we also want to rollback simplex to its original state.
   const unsigned initialSnapshot = simplex.getSnapshot();
 
   // Automatically restore the original state when we return.
   auto restoreState = [&]() {
-    b.removeIdRange(FlatAffineConstraints::IdKind::Local, bInitNumLocals,
+    b.removeIdRange(IntegerPolyhedron::IdKind::Local, bInitNumLocals,
                     b.getNumLocalIds());
     b.removeInequalityRange(bInitNumIneqs, b.getNumInequalities());
     b.removeEqualityRange(bInitNumEqs, b.getNumEqualities());
     simplex.rollback(initialSnapshot);
   };
 
   // Find out which inequalities of sI correspond to division inequalities for
   // the local variables of sI.
for (Optional<std::pair<unsigned, unsigned>> &maybePair : repr) {
     isDivInequality[maybePair->second] = true;
   }
 
   unsigned offset = simplex.getNumConstraints();
   unsigned numLocalsAdded = b.getNumLocalIds() - bInitNumLocals;
   simplex.appendVariable(numLocalsAdded);
 
   unsigned snapshotBeforeIntersect = simplex.getSnapshot();
-  simplex.intersectFlatAffineConstraints(sI);
+  simplex.intersectIntegerPolyhedron(sI);
 
   if (simplex.isEmpty()) {
     /// b ^ s_i is empty, so b \ s_i = b. We move directly to i + 1.
     /// We are ignoring level i completely, so we restore the state
     /// *before* going to level i + 1.
     restoreState();
     subtractRecursively(b, simplex, s, i + 1, result);
     return;
static void subtractRecursively(IntegerPolyhedron &b, Simplex &simplex,
   restoreState();
 }
 
 /// Return the set difference fac \ set.
 ///
 /// The FAC here is modified in subtractRecursively, so it cannot be a const
 /// reference even though it is restored to its original state before returning
 /// from that function.
-PresburgerSet PresburgerSet::getSetDifference(FlatAffineConstraints fac,
+PresburgerSet PresburgerSet::getSetDifference(IntegerPolyhedron fac,
                                               const PresburgerSet &set) {
   assertDimensionsCompatible(fac, set);
   if (fac.isEmptyByGCDTest())
     return PresburgerSet::getEmptySet(fac.getNumDimIds(),
                                       fac.getNumSymbolIds());
 
   PresburgerSet result(fac.getNumDimIds(), fac.getNumSymbolIds());
   Simplex simplex(fac);
   subtractRecursively(fac, simplex, set, 0, result);
   return result;
 }
 
 /// Return the complement of this set.
 PresburgerSet PresburgerSet::complement() const {
   return getSetDifference(
-      FlatAffineConstraints::getUniverse(getNumDims(), getNumSyms()), *this);
+      IntegerPolyhedron::getUniverse(getNumDims(), getNumSyms()), *this);
 }
 
 /// Return the result of subtract the given set from this set, i.e.,
 /// return `this \ set`.
 PresburgerSet PresburgerSet::subtract(const PresburgerSet &set) const {
   assertDimensionsCompatible(set, *this);
   PresburgerSet result(nDim, nSym);
   // We compute (U_i t_i) \ (U_i set_i) as U_i (t_i \ V_i set_i).
-  for (const FlatAffineConstraints &fac : flatAffineConstraints)
+  for (const IntegerPolyhedron &fac : flatAffineConstraints)
     result.unionSetInPlace(getSetDifference(fac, set));
   return result;
 }
 
 /// Two sets S and T are equal iff S contains T and T contains S.
 /// By "S contains T", we mean that S is a superset of or equal to T.
 ///
 /// S contains T iff T \ S is empty, since if T \ S contains a
 /// point then this is a point that is contained in T but not S.
 ///
 /// Therefore, S is equal to T iff S \ T and T \ S are both empty.
 bool PresburgerSet::isEqual(const PresburgerSet &set) const {
   assertDimensionsCompatible(set, *this);
   return this->subtract(set).isIntegerEmpty() &&
          set.subtract(*this).isIntegerEmpty();
 }
 
 /// Return true if all the sets in the union are known to be integer empty,
 /// false otherwise.
 bool PresburgerSet::isIntegerEmpty() const {
   // The set is empty iff all of the disjuncts are empty.
-  for (const FlatAffineConstraints &fac : flatAffineConstraints) {
+  for (const IntegerPolyhedron &fac : flatAffineConstraints) {
     if (!fac.isIntegerEmpty())
       return false;
   }
   return true;
 }
 
 bool PresburgerSet::findIntegerSample(SmallVectorImpl<int64_t> &sample) {
   // A sample exists iff any of the disjuncts contains a sample.
-  for (const FlatAffineConstraints &fac : flatAffineConstraints) {
+  for (const IntegerPolyhedron &fac : flatAffineConstraints) {
     if (Optional<SmallVector<int64_t, 8>> opt = fac.findIntegerSample()) {
       sample = std::move(*opt);
       return true;
     }
   }
   return false;
 }
 
 PresburgerSet PresburgerSet::coalesce() const {
   PresburgerSet newSet = PresburgerSet::getEmptySet(getNumDims(), getNumSyms());
   llvm::SmallBitVector isRedundant(getNumFACs());
 
   for (unsigned i = 0, e = flatAffineConstraints.size(); i < e; ++i) {
     if (isRedundant[i])
       continue;
     Simplex simplex(flatAffineConstraints[i]);
 
     // Check whether the polytope of `simplex` is empty. If so, it is trivially
     // redundant.
     if (simplex.isEmpty()) {
       isRedundant[i] = true;
       continue;
     }
 
-    // Check whether `FlatAffineConstraints[i]` is contained in any FAC, that is
+    // Check whether `IntegerPolyhedron[i]` is contained in any FAC, that is
     // different from itself and not yet marked as redundant.
     for (unsigned j = 0, e = flatAffineConstraints.size(); j < e; ++j) {
       if (j == i || isRedundant[j])
         continue;
 
       if (simplex.isRationalSubsetOf(flatAffineConstraints[j])) {
         isRedundant[i] = true;
         break;
       }
     }
   }
 
   for (unsigned i = 0, e = flatAffineConstraints.size(); i < e; ++i)
     if (!isRedundant[i])
       newSet.unionFACInPlace(flatAffineConstraints[i]);
 
   return newSet;
 }
 
 void PresburgerSet::print(raw_ostream &os) const {
-  os << getNumFACs() << " FlatAffineConstraints:\n";
-  for (const FlatAffineConstraints &fac : flatAffineConstraints) {
+  os << getNumFACs() << " IntegerPolyhedron:\n";
+  for (const IntegerPolyhedron &fac : flatAffineConstraints) {
     fac.print(os);
     os << '\n';
   }
 }
 
 void PresburgerSet::dump() const { print(llvm::errs()); }

mlir/unittests/Analysis/CMakeLists.txt

 add_mlir_unittest(MLIRAnalysisTests
   AffineStructuresParser.cpp
   AffineStructuresParserTest.cpp
   AffineStructuresTest.cpp
-  LinearTransformTest.cpp
   PresburgerSetTest.cpp
 )
 
 target_link_libraries(MLIRAnalysisTests
   PRIVATE
   MLIRLoopAnalysis
   MLIRParser
   )
 
 add_subdirectory(Presburger)

mlir/unittests/Analysis/Presburger/CMakeLists.txt

 add_mlir_unittest(MLIRPresburgerTests
   IntegerPolyhedronTest.cpp
+  LinearTransformTest.cpp
   MatrixTest.cpp
   SimplexTest.cpp
   ../AffineStructuresParser.cpp
 )
 
 target_link_libraries(MLIRPresburgerTests
   PRIVATE MLIRPresburger
   MLIRLoopAnalysis
   MLIRParser
   )

mlir/unittests/Analysis/Presburger/LinearTransformTest.cpp

 //===- LinearTransformTest.cpp - Tests for LinearTransform ----------------===//
 //
 // 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/Presburger/LinearTransform.h"
 #include <gmock/gmock.h>
 #include <gtest/gtest.h>
 
 namespace mlir {
 
 void testColumnEchelonForm(const Matrix &m, unsigned expectedRank) {
   unsigned lastAllowedNonZeroCol = 0;
   std::pair<unsigned, LinearTransform> result =