diff --git a/mlir/include/mlir/Analysis/AffineStructures.h b/mlir/include/mlir/Analysis/AffineStructures.h --- a/mlir/include/mlir/Analysis/AffineStructures.h +++ b/mlir/include/mlir/Analysis/AffineStructures.h @@ -102,65 +102,14 @@ return cst->getKind() == Kind::FlatAffineConstraints; } - /// 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> findIntegerSample() const; - - /// Returns true if the given point satisfies the constraints, or false - /// otherwise. - bool containsPoint(ArrayRef 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> ÷nds, - SmallVector &denominators, - std::vector>> &repr) const; - void getLocalReprs( - std::vector>> &repr) const; - void getLocalReprs(std::vector> ÷nds, - SmallVector &denominators) const; - // Clones this object. std::unique_ptr clone() const; - /// The type of bound: equal, lower bound or upper bound. - enum BoundType { EQ, LB, UB }; + /// The inherited base class functions are not candidates for overload + /// resolution when another function with the same name exists in the derived + /// class (function name hiding), so we have to specify usage of `addBound` + /// here. + using IntegerPolyhedron::addBound; /// 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 @@ -171,12 +120,6 @@ /// 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 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 @@ -193,91 +136,15 @@ SmallVectorImpl *lbMaps, SmallVectorImpl *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 dividend, int64_t divisor); - /// 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); } - - /// 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`. void clearAndCopyFrom(const IntegerPolyhedron &other) override; - /// 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 getConstantBoundOnDimSize( - unsigned pos, SmallVectorImpl *lb = nullptr, - int64_t *boundFloorDivisor = nullptr, - SmallVectorImpl *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 getConstantBound(BoundType type, unsigned pos) const; - /// 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 @@ -290,61 +157,7 @@ unsigned symStartPos, ArrayRef 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 >= 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); - protected: - /// 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 - Optional 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 @@ -358,73 +171,10 @@ LogicalResult flattenAlignedMapAndMergeLocals( AffineMap map, std::vector> *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; - /// Prints the number of constraints, dimensions, symbols and locals in the /// FlatAffineConstraints. Also, prints for each identifier whether there is /// an SSA Value attached to it. void printSpace(raw_ostream &os) const override; - - /// 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; }; /// An extension of FlatAffineConstraints in which dimensions and symbols can diff --git a/mlir/include/mlir/Analysis/Presburger/IntegerPolyhedron.h b/mlir/include/mlir/Analysis/Presburger/IntegerPolyhedron.h --- a/mlir/include/mlir/Analysis/Presburger/IntegerPolyhedron.h +++ b/mlir/include/mlir/Analysis/Presburger/IntegerPolyhedron.h @@ -14,6 +14,7 @@ #define MLIR_ANALYSIS_PRESBURGER_INTEGERPOLYHEDRON_H #include "mlir/Analysis/Presburger/Matrix.h" +#include "mlir/Support/LogicalResult.h" namespace mlir { @@ -214,10 +215,254 @@ SmallVectorImpl *eqIndices = nullptr, unsigned offset = 0, unsigned num = 0) const; + /// 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> findIntegerSample() const; + + /// Returns true if the given point satisfies the constraints, or false + /// otherwise. + bool containsPoint(ArrayRef 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> ÷nds, + SmallVector &denominators, + std::vector>> &repr) const; + void getLocalReprs( + std::vector>> &repr) const; + void getLocalReprs(std::vector> ÷nds, + SmallVector &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 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 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); } + + /// 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 getConstantBoundOnDimSize( + unsigned pos, SmallVectorImpl *lb = nullptr, + int64_t *boundFloorDivisor = nullptr, + SmallVectorImpl *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 getConstantBound(BoundType type, unsigned pos) 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 >= 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: + /// 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 + Optional 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; + /// 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. @@ -238,6 +483,18 @@ /// arrays as needed. virtual void removeIdRange(unsigned idStart, unsigned idLimit); + /// 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; + /// Total number of identifiers. unsigned numIds; diff --git a/mlir/lib/Analysis/AffineStructures.cpp b/mlir/lib/Analysis/AffineStructures.cpp --- a/mlir/lib/Analysis/AffineStructures.cpp +++ b/mlir/lib/Analysis/AffineStructures.cpp @@ -31,8 +31,6 @@ #define DEBUG_TYPE "affine-structures" using namespace mlir; -using llvm::SmallDenseMap; -using llvm::SmallDenseSet; namespace { @@ -698,627 +696,17 @@ 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 -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(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(this, i); - } - for (unsigned i = 0, e = getNumInequalities(); i < e; ++i) { - normalizeConstraintByGCD(this, i); - } -} - 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 times 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 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> -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 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. - 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> 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.getNumIds(); ++j) { - int64_t coeff = cone.atIneq(i, j); - if (coeff < 0) - cone.atIneq(i, cone.getNumIds()) += 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 shrunkenConeSample = - shrunkenConeSimplex.getRationalSample(); - - SmallVector coneSample(llvm::map_range(shrunkenConeSample, ceil)); - - // 6) Return transform * concat(boundedSample, coneSample). - SmallVector &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 expr, ArrayRef 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 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>> &repr) const { - std::vector> dividends(getNumLocalIds()); - SmallVector denominators(getNumLocalIds()); - getLocalReprs(dividends, denominators, repr); -} - -void FlatAffineConstraints::getLocalReprs( - std::vector> ÷nds, - SmallVector &denominators) const { - std::vector>> repr( - getNumLocalIds()); - getLocalReprs(dividends, denominators, repr); -} - -void FlatAffineConstraints::getLocalReprs( - std::vector> ÷nds, - SmallVector &denominators, - std::vector>> &repr) const { - - repr.resize(getNumLocalIds()); - dividends.resize(getNumLocalIds()); - denominators.resize(getNumLocalIds()); - - SmallVector 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(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(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(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: @@ -1481,278 +869,6 @@ 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 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> divsA, divsB; - SmallVector 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 &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(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(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 FlatAffineConstraints::getLowerAndUpperBound( unsigned pos, unsigned offset, unsigned num, unsigned symStartPos, ArrayRef localExprs, MLIRContext *context) const { @@ -2179,61 +1295,6 @@ 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 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 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 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()); - bound[bound.size() - 1] += divisor - 1; - addInequality(bound); -} - bool FlatAffineValueConstraints::findId(Value val, unsigned *pos) const { unsigned i = 0; for (const auto &mayBeId : values) { @@ -2257,13 +1318,6 @@ 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; @@ -2273,289 +1327,6 @@ 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; -} - -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 FlatAffineConstraints::getConstantBoundOnDimSize( - unsigned pos, SmallVectorImpl *lb, int64_t *boundFloorDivisor, - SmallVectorImpl *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 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 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(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 -Optional -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 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 FlatAffineConstraints::getConstantBound(BoundType type, - unsigned pos) const { - assert(type != BoundType::EQ && "EQ not implemented"); - FlatAffineConstraints tmpCst(*this); - if (type == BoundType::LB) - return tmpCst.computeConstantLowerOrUpperBound(pos); - return tmpCst.computeConstantLowerOrUpperBound(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::printSpace(raw_ostream &os) const { IntegerPolyhedron::printSpace(os); os << "("; @@ -2572,77 +1343,6 @@ os << " const)\n"; } -/// 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 >= 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 - // for a given row. - SmallDenseMap, std::pair> - rowsWithoutConstTerm; - // To unique rows. - SmallDenseSet, 8> rowSet; - - // Check if constraint is of the form >= 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 redunIneq(getNumInequalities(), false); - for (unsigned r = 0, e = getNumInequalities(); r < e; r++) { - int64_t *rowStart = &inequalities(r, 0); - auto row = ArrayRef(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(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 IntegerPolyhedron &other) { if (auto *otherValueSet = dyn_cast(&other)) assert(!otherValueSet->hasValues() && @@ -2674,232 +1374,6 @@ values.resize(numIds, None); } -static std::pair -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 lbIndices; - // Positions of constraints that are lower bounds on the variable. - SmallVector ubIndices; - // Positions of constraints that do not involve the variable. - std::vector 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 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 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 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, 8> newVals; @@ -2913,41 +1387,6 @@ 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); @@ -2956,167 +1395,6 @@ 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 a, ArrayRef 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> boundingLbs; - std::vector> boundingUbs; - boundingLbs.reserve(2 * getNumDimIds()); - boundingUbs.reserve(2 * getNumDimIds()); - - // To hold lower and upper bounds for each dimension. - SmallVector lb, otherLb, ub, otherUb; - // To compute min of lower bounds and max of upper bounds for each dimension. - SmallVector minLb(getNumSymbolIds() + 1); - SmallVector maxUb(getNumSymbolIds() + 1); - // To compute final new lower and upper bounds for the union. - SmallVector 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()); - 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"); @@ -3215,12 +1493,6 @@ 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. @@ -3274,55 +1546,6 @@ 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 &nbIneqIndices, - SmallVectorImpl &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 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 *newSyms) { diff --git a/mlir/lib/Analysis/Presburger/IntegerPolyhedron.cpp b/mlir/lib/Analysis/Presburger/IntegerPolyhedron.cpp --- a/mlir/lib/Analysis/Presburger/IntegerPolyhedron.cpp +++ b/mlir/lib/Analysis/Presburger/IntegerPolyhedron.cpp @@ -11,8 +11,18 @@ //===----------------------------------------------------------------------===// #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 "llvm/ADT/DenseMap.h" +#include "llvm/ADT/DenseSet.h" +#include "llvm/Support/Debug.h" + +#define DEBUG_TYPE "presburger" using namespace mlir; +using llvm::SmallDenseMap; +using llvm::SmallDenseSet; std::unique_ptr IntegerPolyhedron::clone() const { return std::make_unique(*this); @@ -307,6 +317,1778 @@ *this = other; } +// 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 +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(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(this, i); + } + for (unsigned i = 0, e = getNumInequalities(); i < e; ++i) { + normalizeConstraintByGCD(this, i); + } +} + +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 times 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 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 &poly, unsigned eqIndex, + unsigned numDims) { + for (unsigned e = poly.getNumIds(), j = e - numDims; j < e; ++j) + if (poly.atEq(eqIndex, j) != 0) + return true; + return false; +} +bool ineqInvolvesSuffixDims(const IntegerPolyhedron &poly, unsigned ineqIndex, + unsigned numDims) { + for (unsigned e = poly.getNumIds(), j = e - numDims; j < e; ++j) + if (poly.atIneq(ineqIndex, j) != 0) + return true; + return false; +} + +void removeConstraintsInvolvingSuffixDims(IntegerPolyhedron &poly, + 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 = poly.getNumEqualities(); i > 0; i--) + if (eqInvolvesSuffixDims(poly, i - 1, unboundedDims)) + poly.removeEquality(i - 1); + for (unsigned i = poly.getNumInequalities(); i > 0; i--) + if (ineqInvolvesSuffixDims(poly, i - 1, unboundedDims)) + poly.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> 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 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> 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.getNumIds(); ++j) { + int64_t coeff = cone.atIneq(i, j); + if (coeff < 0) + cone.atIneq(i, cone.getNumIds()) += 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 shrunkenConeSample = + shrunkenConeSimplex.getRationalSample(); + + SmallVector coneSample(llvm::map_range(shrunkenConeSample, ceil)); + + // 6) Return transform * concat(boundedSample, coneSample). + SmallVector &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 expr, ArrayRef 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 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>> &repr) const { + std::vector> dividends(getNumLocalIds()); + SmallVector denominators(getNumLocalIds()); + getLocalReprs(dividends, denominators, repr); +} + +void IntegerPolyhedron::getLocalReprs( + std::vector> ÷nds, + SmallVector &denominators) const { + std::vector>> repr( + getNumLocalIds()); + getLocalReprs(dividends, denominators, repr); +} + +void IntegerPolyhedron::getLocalReprs( + std::vector> ÷nds, + SmallVector &denominators, + std::vector>> &repr) const { + + repr.resize(getNumLocalIds()); + dividends.resize(getNumLocalIds()); + denominators.resize(getNumLocalIds()); + + SmallVector 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(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(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(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 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 &poly, unsigned pos1, + unsigned pos2) { + + assert(pos1 < poly.getNumLocalIds() && "Invalid local id position"); + assert(pos2 < poly.getNumLocalIds() && "Invalid local id position"); + + unsigned localOffset = poly.getNumDimAndSymbolIds(); + pos1 += localOffset; + pos2 += localOffset; + for (unsigned i = 0, e = poly.getNumInequalities(); i < e; ++i) + poly.atIneq(i, pos1) += poly.atIneq(i, pos2); + for (unsigned i = 0, e = poly.getNumEqualities(); i < e; ++i) + poly.atEq(i, pos1) += poly.atEq(i, pos2); + poly.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 &polyA = *this; + IntegerPolyhedron &polyB = other; + + // Merge local ids of polyA and polyB without using division information, + // i.e. append local ids of `polyB` to `polyA` and insert local ids of `polyA` + // to `polyB` at start of its local ids. + unsigned initLocals = polyA.getNumLocalIds(); + insertLocalId(polyA.getNumLocalIds(), polyB.getNumLocalIds()); + polyB.insertLocalId(0, initLocals); + + // Get division representations from each poly. + std::vector> divsA, divsB; + SmallVector denomsA, denomsB; + polyA.getLocalReprs(divsA, denomsA); + polyB.getLocalReprs(divsB, denomsB); + + // Copy division information for polyB into `divsA` and `denomsA`, so that + // these have the combined division information of both polys. Since newly + // added local variables in polyA and polyB 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(polyA, i, j); + eliminateRedundantLocalId(polyB, i, j); + for (unsigned k = 0, g = divsA.size(); k < g; ++k) { + SmallVector &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(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(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 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 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 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()); + 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; +} + +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 IntegerPolyhedron::getConstantBoundOnDimSize( + unsigned pos, SmallVectorImpl *lb, int64_t *boundFloorDivisor, + SmallVectorImpl *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 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 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(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 +Optional +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 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 IntegerPolyhedron::getConstantBound(BoundType type, + unsigned pos) const { + assert(type != BoundType::EQ && "EQ not implemented"); + IntegerPolyhedron tmpCst(*this); + if (type == BoundType::LB) + return tmpCst.computeConstantLowerOrUpperBound(pos); + return tmpCst.computeConstantLowerOrUpperBound(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 >= 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 + // for a given row. + SmallDenseMap, std::pair> + rowsWithoutConstTerm; + // To unique rows. + SmallDenseSet, 8> rowSet; + + // Check if constraint is of the form >= 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 redunIneq(getNumInequalities(), false); + for (unsigned r = 0, e = getNumInequalities(); r < e; r++) { + int64_t *rowStart = &inequalities(r, 0); + auto row = ArrayRef(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(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 +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}; +} + +#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 lbIndices; + // Positions of constraints that are lower bounds on the variable. + SmallVector ubIndices; + // Positions of constraints that do not involve the variable. + std::vector 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 newPoly( + 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 ineq; + ineq.reserve(newPoly.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. + newPoly.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 ineq; + ineq.reserve(getNumCols() - 1); + for (unsigned l = 0, e = getNumCols(); l < e; l++) { + if (l == pos) + continue; + ineq.push_back(atIneq(nbPos, l)); + } + newPoly.addInequality(ineq); + } + + assert(newPoly.getNumConstraints() == + lbIndices.size() * ubIndices.size() + nbIndices.size()); + + // Copy over the equalities. + for (unsigned r = 0, e = getNumEqualities(); r < e; r++) { + SmallVector eq; + eq.reserve(newPoly.getNumCols()); + for (unsigned l = 0, e = getNumCols(); l < e; l++) { + if (l == pos) + continue; + eq.push_back(atEq(r, l)); + } + newPoly.addEquality(eq); + } + + // GCD tightening and normalization allows detection of more trivially + // redundant constraints. + newPoly.gcdTightenInequalities(); + newPoly.normalizeConstraintsByGCD(); + newPoly.removeTrivialRedundancy(); + clearAndCopyFrom(newPoly); + 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 a, ArrayRef 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> boundingLbs; + std::vector> boundingUbs; + boundingLbs.reserve(2 * getNumDimIds()); + boundingUbs.reserve(2 * getNumDimIds()); + + // To hold lower and upper bounds for each dimension. + SmallVector lb, otherLb, ub, otherUb; + // To compute min of lower bounds and max of upper bounds for each dimension. + SmallVector minLb(getNumSymbolIds() + 1); + SmallVector maxUb(getNumSymbolIds() + 1); + // To compute final new lower and upper bounds for the union. + SmallVector 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()); + 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 &nbIneqIndices, + SmallVectorImpl &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 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() diff --git a/mlir/unittests/Analysis/AffineStructuresTest.cpp b/mlir/unittests/Analysis/AffineStructuresTest.cpp deleted file mode 100644 --- a/mlir/unittests/Analysis/AffineStructuresTest.cpp +++ /dev/null @@ -1,830 +0,0 @@ -//===- AffineStructuresTest.cpp - Tests for AffineStructures ----*- 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 -// -//===----------------------------------------------------------------------===// - -#include "mlir/Analysis/AffineStructures.h" -#include "./AffineStructuresParser.h" -#include "mlir/IR/IntegerSet.h" -#include "mlir/IR/MLIRContext.h" - -#include -#include - -#include - -namespace mlir { - -using testing::ElementsAre; - -enum class TestFunction { Sample, Empty }; - -/// If fn is TestFunction::Sample (default): -/// If hasSample is true, check that findIntegerSample returns a valid sample -/// for the FlatAffineConstraints fac. -/// If hasSample is false, check that findIntegerSample returns None. -/// -/// If fn is TestFunction::Empty, check that isIntegerEmpty returns the -/// opposite of hasSample. -static void checkSample(bool hasSample, const FlatAffineConstraints &fac, - TestFunction fn = TestFunction::Sample) { - Optional> maybeSample; - switch (fn) { - case TestFunction::Sample: - maybeSample = fac.findIntegerSample(); - if (!hasSample) { - EXPECT_FALSE(maybeSample.hasValue()); - if (maybeSample.hasValue()) { - for (auto x : *maybeSample) - llvm::errs() << x << ' '; - llvm::errs() << '\n'; - } - } else { - ASSERT_TRUE(maybeSample.hasValue()); - EXPECT_TRUE(fac.containsPoint(*maybeSample)); - } - break; - case TestFunction::Empty: - EXPECT_EQ(!hasSample, fac.isIntegerEmpty()); - break; - } -} - -/// Construct a FlatAffineConstraints from a set of inequality and -/// equality constraints. -static FlatAffineConstraints -makeFACFromConstraints(unsigned ids, ArrayRef> ineqs, - ArrayRef> eqs, - unsigned syms = 0) { - FlatAffineConstraints fac(ineqs.size(), eqs.size(), ids + 1, ids - syms, syms, - /*numLocals=*/0); - for (const auto &eq : eqs) - fac.addEquality(eq); - for (const auto &ineq : ineqs) - fac.addInequality(ineq); - return fac; -} - -/// Check sampling for all the permutations of the dimensions for the given -/// constraint set. Since the GBR algorithm progresses dimension-wise, different -/// orderings may cause the algorithm to proceed differently. At least some of -///.these permutations should make it past the heuristics and test the -/// implementation of the GBR algorithm itself. -/// Use TestFunction fn to test. -static void checkPermutationsSample(bool hasSample, unsigned nDim, - ArrayRef> ineqs, - ArrayRef> eqs, - TestFunction fn = TestFunction::Sample) { - SmallVector perm(nDim); - std::iota(perm.begin(), perm.end(), 0); - auto permute = [&perm](ArrayRef coeffs) { - SmallVector permuted; - for (unsigned id : perm) - permuted.push_back(coeffs[id]); - permuted.push_back(coeffs.back()); - return permuted; - }; - do { - SmallVector, 4> permutedIneqs, permutedEqs; - for (const auto &ineq : ineqs) - permutedIneqs.push_back(permute(ineq)); - for (const auto &eq : eqs) - permutedEqs.push_back(permute(eq)); - - checkSample(hasSample, - makeFACFromConstraints(nDim, permutedIneqs, permutedEqs), fn); - } while (std::next_permutation(perm.begin(), perm.end())); -} - -/// Parses a FlatAffineConstraints from a StringRef. It is expected that the -/// string represents a valid IntegerSet, otherwise it will violate a gtest -/// assertion. -static FlatAffineConstraints parseFAC(StringRef str, MLIRContext *context) { - FailureOr fac = parseIntegerSetToFAC(str, context); - - EXPECT_TRUE(succeeded(fac)); - - return *fac; -} - -TEST(FlatAffineConstraintsTest, FindSampleTest) { - // Bounded sets with only inequalities. - - MLIRContext context; - - // 0 <= 7x <= 5 - checkSample(true, parseFAC("(x) : (7 * x >= 0, -7 * x + 5 >= 0)", &context)); - - // 1 <= 5x and 5x <= 4 (no solution). - checkSample(false, - parseFAC("(x) : (5 * x - 1 >= 0, -5 * x + 4 >= 0)", &context)); - - // 1 <= 5x and 5x <= 9 (solution: x = 1). - checkSample(true, - parseFAC("(x) : (5 * x - 1 >= 0, -5 * x + 9 >= 0)", &context)); - - // Bounded sets with equalities. - // x >= 8 and 40 >= y and x = y. - checkSample(true, parseFAC("(x,y) : (x - 8 >= 0, -y + 40 >= 0, x - y == 0)", - &context)); - - // x <= 10 and y <= 10 and 10 <= z and x + 2y = 3z. - // solution: x = y = z = 10. - checkSample(true, parseFAC("(x,y,z) : (-x + 10 >= 0, -y + 10 >= 0, " - "z - 10 >= 0, x + 2 * y - 3 * z == 0)", - &context)); - - // x <= 10 and y <= 10 and 11 <= z and x + 2y = 3z. - // This implies x + 2y >= 33 and x + 2y <= 30, which has no solution. - checkSample(false, parseFAC("(x,y,z) : (-x + 10 >= 0, -y + 10 >= 0, " - "z - 11 >= 0, x + 2 * y - 3 * z == 0)", - &context)); - - // 0 <= r and r <= 3 and 4q + r = 7. - // Solution: q = 1, r = 3. - checkSample( - true, - parseFAC("(q,r) : (r >= 0, -r + 3 >= 0, 4 * q + r - 7 == 0)", &context)); - - // 4q + r = 7 and r = 0. - // Solution: q = 1, r = 3. - checkSample(false, - parseFAC("(q,r) : (4 * q + r - 7 == 0, r == 0)", &context)); - - // The next two sets are large sets that should take a long time to sample - // with a naive branch and bound algorithm but can be sampled efficiently with - // the GBR algorithm. - // - // This is a triangle with vertices at (1/3, 0), (2/3, 0) and (10000, 10000). - checkSample(true, parseFAC("(x,y) : (y >= 0, " - "300000 * x - 299999 * y - 100000 >= 0, " - "-300000 * x + 299998 * y + 200000 >= 0)", - &context)); - - // This is a tetrahedron with vertices at - // (1/3, 0, 0), (2/3, 0, 0), (2/3, 0, 10000), and (10000, 10000, 10000). - // The first three points form a triangular base on the xz plane with the - // apex at the fourth point, which is the only integer point. - checkPermutationsSample( - true, 3, - { - {0, 1, 0, 0}, // y >= 0 - {0, -1, 1, 0}, // z >= y - {300000, -299998, -1, - -100000}, // -300000x + 299998y + 100000 + z <= 0. - {-150000, 149999, 0, 100000}, // -150000x + 149999y + 100000 >= 0. - }, - {}); - - // Same thing with some spurious extra dimensions equated to constants. - checkSample( - true, - parseFAC("(a,b,c,d,e) : (b + d - e >= 0, -b + c - d + e >= 0, " - "300000 * a - 299998 * b - c - 9 * d + 21 * e - 112000 >= 0, " - "-150000 * a + 149999 * b - 15 * d + 47 * e + 68000 >= 0, " - "d - e == 0, d + e - 2000 == 0)", - &context)); - - // This is a tetrahedron with vertices at - // (1/3, 0, 0), (2/3, 0, 0), (2/3, 0, 100), (100, 100 - 1/3, 100). - checkPermutationsSample(false, 3, - { - {0, 1, 0, 0}, - {0, -300, 299, 0}, - {300 * 299, -89400, -299, -100 * 299}, - {-897, 894, 0, 598}, - }, - {}); - - // Two tests involving equalities that are integer empty but not rational - // empty. - - // This is a line segment from (0, 1/3) to (100, 100 + 1/3). - checkSample( - false, parseFAC("(x,y) : (x >= 0, -x + 100 >= 0, 3 * x - 3 * y + 1 == 0)", - &context)); - - // A thin parallelogram. 0 <= x <= 100 and x + 1/3 <= y <= x + 2/3. - checkSample(false, - parseFAC("(x,y) : (x >= 0, -x + 100 >= 0, " - "3 * x - 3 * y + 2 >= 0, -3 * x + 3 * y - 1 >= 0)", - &context)); - - checkSample(true, parseFAC("(x,y) : (2 * x >= 0, -2 * x + 99 >= 0, " - "2 * y >= 0, -2 * y + 99 >= 0)", - &context)); - - // 2D cone with apex at (10000, 10000) and - // edges passing through (1/3, 0) and (2/3, 0). - checkSample(true, parseFAC("(x,y) : (300000 * x - 299999 * y - 100000 >= 0, " - "-300000 * x + 299998 * y + 200000 >= 0)", - &context)); - - // Cartesian product of a tetrahedron and a 2D cone. - // The tetrahedron has vertices at - // (1/3, 0, 0), (2/3, 0, 0), (2/3, 0, 10000), and (10000, 10000, 10000). - // The first three points form a triangular base on the xz plane with the - // apex at the fourth point, which is the only integer point. - // The cone has apex at (10000, 10000) and - // edges passing through (1/3, 0) and (2/3, 0). - checkPermutationsSample( - true /* not empty */, 5, - { - // Tetrahedron contraints: - {0, 1, 0, 0, 0, 0}, // y >= 0 - {0, -1, 1, 0, 0, 0}, // z >= y - // -300000x + 299998y + 100000 + z <= 0. - {300000, -299998, -1, 0, 0, -100000}, - // -150000x + 149999y + 100000 >= 0. - {-150000, 149999, 0, 0, 0, 100000}, - - // Triangle constraints: - // 300000p - 299999q >= 100000 - {0, 0, 0, 300000, -299999, -100000}, - // -300000p + 299998q + 200000 >= 0 - {0, 0, 0, -300000, 299998, 200000}, - }, - {}); - - // Cartesian product of same tetrahedron as above and {(p, q) : 1/3 <= p <= - // 2/3}. Since the second set is empty, the whole set is too. - checkPermutationsSample( - false /* empty */, 5, - { - // Tetrahedron contraints: - {0, 1, 0, 0, 0, 0}, // y >= 0 - {0, -1, 1, 0, 0, 0}, // z >= y - // -300000x + 299998y + 100000 + z <= 0. - {300000, -299998, -1, 0, 0, -100000}, - // -150000x + 149999y + 100000 >= 0. - {-150000, 149999, 0, 0, 0, 100000}, - - // Second set constraints: - // 3p >= 1 - {0, 0, 0, 3, 0, -1}, - // 3p <= 2 - {0, 0, 0, -3, 0, 2}, - }, - {}); - - // Cartesian product of same tetrahedron as above and - // {(p, q, r) : 1 <= p <= 2 and p = 3q + 3r}. - // Since the second set is empty, the whole set is too. - checkPermutationsSample( - false /* empty */, 5, - { - // Tetrahedron contraints: - {0, 1, 0, 0, 0, 0, 0}, // y >= 0 - {0, -1, 1, 0, 0, 0, 0}, // z >= y - // -300000x + 299998y + 100000 + z <= 0. - {300000, -299998, -1, 0, 0, 0, -100000}, - // -150000x + 149999y + 100000 >= 0. - {-150000, 149999, 0, 0, 0, 0, 100000}, - - // Second set constraints: - // p >= 1 - {0, 0, 0, 1, 0, 0, -1}, - // p <= 2 - {0, 0, 0, -1, 0, 0, 2}, - }, - { - {0, 0, 0, 1, -3, -3, 0}, // p = 3q + 3r - }); - - // Cartesian product of a tetrahedron and a 2D cone. - // The tetrahedron is empty and has vertices at - // (1/3, 0, 0), (2/3, 0, 0), (2/3, 0, 100), and (100, 100 - 1/3, 100). - // The cone has apex at (10000, 10000) and - // edges passing through (1/3, 0) and (2/3, 0). - // Since the tetrahedron is empty, the Cartesian product is too. - checkPermutationsSample(false /* empty */, 5, - { - // Tetrahedron contraints: - {0, 1, 0, 0, 0, 0}, - {0, -300, 299, 0, 0, 0}, - {300 * 299, -89400, -299, 0, 0, -100 * 299}, - {-897, 894, 0, 0, 0, 598}, - - // Triangle constraints: - // 300000p - 299999q >= 100000 - {0, 0, 0, 300000, -299999, -100000}, - // -300000p + 299998q + 200000 >= 0 - {0, 0, 0, -300000, 299998, 200000}, - }, - {}); - - // Cartesian product of same tetrahedron as above and - // {(p, q) : 1/3 <= p <= 2/3}. - checkPermutationsSample(false /* empty */, 5, - { - // Tetrahedron contraints: - {0, 1, 0, 0, 0, 0}, - {0, -300, 299, 0, 0, 0}, - {300 * 299, -89400, -299, 0, 0, -100 * 299}, - {-897, 894, 0, 0, 0, 598}, - - // Second set constraints: - // 3p >= 1 - {0, 0, 0, 3, 0, -1}, - // 3p <= 2 - {0, 0, 0, -3, 0, 2}, - }, - {}); - - checkSample(true, parseFAC("(x, y, z) : (2 * x - 1 >= 0, x - y - 1 == 0, " - "y - z == 0)", - &context)); - - // Regression tests for the computation of dual coefficients. - checkSample(false, parseFAC("(x, y, z) : (" - "6*x - 4*y + 9*z + 2 >= 0," - "x + 5*y + z + 5 >= 0," - "-4*x + y + 2*z - 1 >= 0," - "-3*x - 2*y - 7*z - 1 >= 0," - "-7*x - 5*y - 9*z - 1 >= 0)", - &context)); - checkSample(true, parseFAC("(x, y, z) : (" - "3*x + 3*y + 3 >= 0," - "-4*x - 8*y - z + 4 >= 0," - "-7*x - 4*y + z + 1 >= 0," - "2*x - 7*y - 8*z - 7 >= 0," - "9*x + 8*y - 9*z - 7 >= 0)", - &context)); -} - -TEST(FlatAffineConstraintsTest, IsIntegerEmptyTest) { - - MLIRContext context; - - // 1 <= 5x and 5x <= 4 (no solution). - EXPECT_TRUE(parseFAC("(x) : (5 * x - 1 >= 0, -5 * x + 4 >= 0)", &context) - .isIntegerEmpty()); - // 1 <= 5x and 5x <= 9 (solution: x = 1). - EXPECT_FALSE(parseFAC("(x) : (5 * x - 1 >= 0, -5 * x + 9 >= 0)", &context) - .isIntegerEmpty()); - - // Unbounded sets. - EXPECT_TRUE(parseFAC("(x,y,z) : (2 * y - 1 >= 0, -2 * y + 1 >= 0, " - "2 * z - 1 >= 0, 2 * x - 1 == 0)", - &context) - .isIntegerEmpty()); - - EXPECT_FALSE(parseFAC("(x,y,z) : (2 * x - 1 >= 0, -3 * x + 3 >= 0, " - "5 * z - 6 >= 0, -7 * z + 17 >= 0, 3 * y - 2 >= 0)", - &context) - .isIntegerEmpty()); - - EXPECT_FALSE( - parseFAC("(x,y,z) : (2 * x - 1 >= 0, x - y - 1 == 0, y - z == 0)", - &context) - .isIntegerEmpty()); - - // FlatAffineConstraints::isEmpty() does not detect the following sets to be - // empty. - - // 3x + 7y = 1 and 0 <= x, y <= 10. - // Since x and y are non-negative, 3x + 7y can never be 1. - EXPECT_TRUE(parseFAC("(x,y) : (x >= 0, -x + 10 >= 0, y >= 0, -y + 10 >= 0, " - "3 * x + 7 * y - 1 == 0)", - &context) - .isIntegerEmpty()); - - // 2x = 3y and y = x - 1 and x + y = 6z + 2 and 0 <= x, y <= 100. - // Substituting y = x - 1 in 3y = 2x, we obtain x = 3 and hence y = 2. - // Since x + y = 5 cannot be equal to 6z + 2 for any z, the set is empty. - EXPECT_TRUE( - parseFAC("(x,y,z) : (x >= 0, -x + 100 >= 0, y >= 0, -y + 100 >= 0, " - "2 * x - 3 * y == 0, x - y - 1 == 0, x + y - 6 * z - 2 == 0)", - &context) - .isIntegerEmpty()); - - // 2x = 3y and y = x - 1 + 6z and x + y = 6q + 2 and 0 <= x, y <= 100. - // 2x = 3y implies x is a multiple of 3 and y is even. - // Now y = x - 1 + 6z implies y = 2 mod 3. In fact, since y is even, we have - // y = 2 mod 6. Then since x = y + 1 + 6z, we have x = 3 mod 6, implying - // x + y = 5 mod 6, which contradicts x + y = 6q + 2, so the set is empty. - EXPECT_TRUE( - parseFAC( - "(x,y,z,q) : (x >= 0, -x + 100 >= 0, y >= 0, -y + 100 >= 0, " - "2 * x - 3 * y == 0, x - y + 6 * z - 1 == 0, x + y - 6 * q - 2 == 0)", - &context) - .isIntegerEmpty()); - - // Set with symbols. - EXPECT_FALSE( - parseFAC("(x)[s] : (x + s >= 0, x - s == 0)", &context).isIntegerEmpty()); -} - -TEST(FlatAffineConstraintsTest, removeRedundantConstraintsTest) { - MLIRContext context; - - FlatAffineConstraints fac = - parseFAC("(x) : (x - 2 >= 0, -x + 2 >= 0, x - 2 == 0)", &context); - fac.removeRedundantConstraints(); - - // Both inequalities are redundant given the equality. Both have been removed. - EXPECT_EQ(fac.getNumInequalities(), 0u); - EXPECT_EQ(fac.getNumEqualities(), 1u); - - FlatAffineConstraints fac2 = - parseFAC("(x,y) : (x - 3 >= 0, y - 2 >= 0, x - y == 0)", &context); - fac2.removeRedundantConstraints(); - - // The second inequality is redundant and should have been removed. The - // remaining inequality should be the first one. - EXPECT_EQ(fac2.getNumInequalities(), 1u); - EXPECT_THAT(fac2.getInequality(0), ElementsAre(1, 0, -3)); - EXPECT_EQ(fac2.getNumEqualities(), 1u); - - FlatAffineConstraints fac3 = - parseFAC("(x,y,z) : (x - y == 0, x - z == 0, y - z == 0)", &context); - fac3.removeRedundantConstraints(); - - // One of the three equalities can be removed. - EXPECT_EQ(fac3.getNumInequalities(), 0u); - EXPECT_EQ(fac3.getNumEqualities(), 2u); - - FlatAffineConstraints fac4 = - parseFAC("(a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q) : (" - "b - 1 >= 0," - "-b + 500 >= 0," - "-16 * d + f >= 0," - "f - 1 >= 0," - "-f + 998 >= 0," - "16 * d - f + 15 >= 0," - "-16 * e + g >= 0," - "g - 1 >= 0," - "-g + 998 >= 0," - "16 * e - g + 15 >= 0," - "h >= 0," - "-h + 1 >= 0," - "j - 1 >= 0," - "-j + 500 >= 0," - "-f + 16 * l + 15 >= 0," - "f - 16 * l >= 0," - "-16 * m + o >= 0," - "o - 1 >= 0," - "-o + 998 >= 0," - "16 * m - o + 15 >= 0," - "p >= 0," - "-p + 1 >= 0," - "-g - h + 8 * q + 8 >= 0," - "-o - p + 8 * q + 8 >= 0," - "o + p - 8 * q - 1 >= 0," - "g + h - 8 * q - 1 >= 0," - "-f + n >= 0," - "f - n >= 0," - "k - 10 >= 0," - "-k + 10 >= 0," - "i - 13 >= 0," - "-i + 13 >= 0," - "c - 10 >= 0," - "-c + 10 >= 0," - "a - 13 >= 0," - "-a + 13 >= 0" - ")", - &context); - - // The above is a large set of constraints without any redundant constraints, - // as verified by the Fourier-Motzkin based removeRedundantInequalities. - unsigned nIneq = fac4.getNumInequalities(); - unsigned nEq = fac4.getNumEqualities(); - fac4.removeRedundantInequalities(); - ASSERT_EQ(fac4.getNumInequalities(), nIneq); - ASSERT_EQ(fac4.getNumEqualities(), nEq); - // Now we test that removeRedundantConstraints does not find any constraints - // to be redundant either. - fac4.removeRedundantConstraints(); - EXPECT_EQ(fac4.getNumInequalities(), nIneq); - EXPECT_EQ(fac4.getNumEqualities(), nEq); - - FlatAffineConstraints fac5 = parseFAC( - "(x,y) : (128 * x + 127 >= 0, -x + 7 >= 0, -128 * x + y >= 0, y >= 0)", - &context); - // 128x + 127 >= 0 implies that 128x >= 0, since x has to be an integer. - // (This should be caught by GCDTightenInqualities().) - // So -128x + y >= 0 and 128x + 127 >= 0 imply y >= 0 since we have - // y >= 128x >= 0. - fac5.removeRedundantConstraints(); - EXPECT_EQ(fac5.getNumInequalities(), 3u); - SmallVector redundantConstraint = {0, 1, 0}; - for (unsigned i = 0; i < 3; ++i) { - // Ensure that the removed constraint was the redundant constraint [3]. - EXPECT_NE(fac5.getInequality(i), ArrayRef(redundantConstraint)); - } -} - -TEST(FlatAffineConstraintsTest, addConstantUpperBound) { - FlatAffineConstraints fac(2); - fac.addBound(FlatAffineConstraints::UB, 0, 1); - EXPECT_EQ(fac.atIneq(0, 0), -1); - EXPECT_EQ(fac.atIneq(0, 1), 0); - EXPECT_EQ(fac.atIneq(0, 2), 1); - - fac.addBound(FlatAffineConstraints::UB, {1, 2, 3}, 1); - EXPECT_EQ(fac.atIneq(1, 0), -1); - EXPECT_EQ(fac.atIneq(1, 1), -2); - EXPECT_EQ(fac.atIneq(1, 2), -2); -} - -TEST(FlatAffineConstraintsTest, addConstantLowerBound) { - FlatAffineConstraints fac(2); - fac.addBound(FlatAffineConstraints::LB, 0, 1); - EXPECT_EQ(fac.atIneq(0, 0), 1); - EXPECT_EQ(fac.atIneq(0, 1), 0); - EXPECT_EQ(fac.atIneq(0, 2), -1); - - fac.addBound(FlatAffineConstraints::LB, {1, 2, 3}, 1); - EXPECT_EQ(fac.atIneq(1, 0), 1); - EXPECT_EQ(fac.atIneq(1, 1), 2); - EXPECT_EQ(fac.atIneq(1, 2), 2); -} - -/// Check if the expected division representation of local variables matches the -/// computed representation. The expected division representation is given as -/// a vector of expressions set in `expectedDividends` and the corressponding -/// denominator in `expectedDenominators`. The `denominators` and `dividends` -/// obtained through `getLocalRepr` function is verified against the -/// `expectedDenominators` and `expectedDividends` respectively. -static void checkDivisionRepresentation( - FlatAffineConstraints &fac, - const std::vector> &expectedDividends, - const SmallVectorImpl &expectedDenominators) { - - std::vector> dividends; - SmallVector denominators; - - fac.getLocalReprs(dividends, denominators); - - // Check that the `denominators` and `expectedDenominators` match. - EXPECT_TRUE(expectedDenominators == denominators); - - // Check that the `dividends` and `expectedDividends` match. If the - // denominator for a division is zero, we ignore its dividend. - EXPECT_TRUE(dividends.size() == expectedDividends.size()); - for (unsigned i = 0, e = dividends.size(); i < e; ++i) - if (denominators[i] != 0) - EXPECT_TRUE(expectedDividends[i] == dividends[i]); -} - -TEST(FlatAffineConstraintsTest, computeLocalReprSimple) { - FlatAffineConstraints fac(1); - - fac.addLocalFloorDiv({1, 4}, 10); - fac.addLocalFloorDiv({1, 0, 100}, 10); - - std::vector> divisions = {{1, 0, 0, 4}, - {1, 0, 0, 100}}; - SmallVector denoms = {10, 10}; - - // Check if floordivs can be computed when no other inequalities exist - // and floor divs do not depend on each other. - checkDivisionRepresentation(fac, divisions, denoms); -} - -TEST(FlatAffineConstraintsTest, computeLocalReprConstantFloorDiv) { - FlatAffineConstraints fac(4); - - fac.addInequality({1, 0, 3, 1, 2}); - fac.addInequality({1, 2, -8, 1, 10}); - fac.addEquality({1, 2, -4, 1, 10}); - - fac.addLocalFloorDiv({0, 0, 0, 0, 10}, 30); - fac.addLocalFloorDiv({0, 0, 0, 0, 0, 99}, 101); - - std::vector> divisions = {{0, 0, 0, 0, 0, 0, 10}, - {0, 0, 0, 0, 0, 0, 99}}; - SmallVector denoms = {30, 101}; - - // Check if floordivs with constant numerator can be computed. - checkDivisionRepresentation(fac, divisions, denoms); -} - -TEST(FlatAffineConstraintsTest, computeLocalReprRecursive) { - FlatAffineConstraints fac(4); - fac.addInequality({1, 0, 3, 1, 2}); - fac.addInequality({1, 2, -8, 1, 10}); - fac.addEquality({1, 2, -4, 1, 10}); - - fac.addLocalFloorDiv({0, -2, 7, 2, 10}, 3); - fac.addLocalFloorDiv({3, 0, 9, 2, 2, 10}, 5); - fac.addLocalFloorDiv({0, 1, -123, 2, 0, -4, 10}, 3); - - fac.addInequality({1, 2, -2, 1, -5, 0, 6, 100}); - fac.addInequality({1, 2, -8, 1, 3, 7, 0, -9}); - - std::vector> divisions = { - {0, -2, 7, 2, 0, 0, 0, 10}, - {3, 0, 9, 2, 2, 0, 0, 10}, - {0, 1, -123, 2, 0, -4, 0, 10}}; - - SmallVector denoms = {3, 5, 3}; - - // Check if floordivs which may depend on other floordivs can be computed. - checkDivisionRepresentation(fac, divisions, denoms); -} - -TEST(FlatAffineConstraintsTest, computeLocalReprTightUpperBound) { - MLIRContext context; - - { - FlatAffineConstraints fac = parseFAC("(i) : (i mod 3 - 1 >= 0)", &context); - - // The set formed by the fac is: - // 3q - i + 2 >= 0 <-- Division lower bound - // -3q + i - 1 >= 0 - // -3q + i >= 0 <-- Division upper bound - // We remove redundant constraints to get the set: - // 3q - i + 2 >= 0 <-- Division lower bound - // -3q + i - 1 >= 0 <-- Tighter division upper bound - // thus, making the upper bound tighter. - fac.removeRedundantConstraints(); - - std::vector> divisions = {{1, 0, 0}}; - SmallVector denoms = {3}; - - // Check if the divisions can be computed even with a tighter upper bound. - checkDivisionRepresentation(fac, divisions, denoms); - } - - { - FlatAffineConstraints fac = parseFAC( - "(i, j, q) : (4*q - i - j + 2 >= 0, -4*q + i + j >= 0)", &context); - // Convert `q` to a local variable. - fac.convertDimToLocal(2, 3); - - std::vector> divisions = {{1, 1, 0, 1}}; - SmallVector denoms = {4}; - - // Check if the divisions can be computed even with a tighter upper bound. - checkDivisionRepresentation(fac, divisions, denoms); - } -} - -TEST(FlatAffineConstraintsTest, computeLocalReprNoRepr) { - MLIRContext context; - FlatAffineConstraints fac = - parseFAC("(x, q) : (x - 3 * q >= 0, -x + 3 * q + 3 >= 0)", &context); - // Convert q to a local variable. - fac.convertDimToLocal(1, 2); - - std::vector> divisions = {{0, 0, 0}}; - SmallVector denoms = {0}; - - // Check that no division is computed. - checkDivisionRepresentation(fac, divisions, denoms); -} - -TEST(FlatAffineConstraintsTest, simplifyLocalsTest) { - // (x) : (exists y: 2x + y = 1 and y = 2). - FlatAffineConstraints fac(1, 0, 1); - fac.addEquality({2, 1, -1}); - fac.addEquality({0, 1, -2}); - - EXPECT_TRUE(fac.isEmpty()); - - // (x) : (exists y, z, w: 3x + y = 1 and 2y = z and 3y = w and z = w). - FlatAffineConstraints fac2(1, 0, 3); - fac2.addEquality({3, 1, 0, 0, -1}); - fac2.addEquality({0, 2, -1, 0, 0}); - fac2.addEquality({0, 3, 0, -1, 0}); - fac2.addEquality({0, 0, 1, -1, 0}); - - EXPECT_TRUE(fac2.isEmpty()); - - // (x) : (exists y: x >= y + 1 and 2x + y = 0 and y >= -1). - FlatAffineConstraints fac3(1, 0, 1); - fac3.addInequality({1, -1, -1}); - fac3.addInequality({0, 1, 1}); - fac3.addEquality({2, 1, 0}); - - EXPECT_TRUE(fac3.isEmpty()); -} - -TEST(FlatAffineConstraintsTest, mergeDivisionsSimple) { - { - // (x) : (exists z, y = [x / 2] : x = 3y and x + z + 1 >= 0). - FlatAffineConstraints fac1(1, 0, 1); - fac1.addLocalFloorDiv({1, 0, 0}, 2); // y = [x / 2]. - fac1.addEquality({1, 0, -3, 0}); // x = 3y. - fac1.addInequality({1, 1, 0, 1}); // x + z + 1 >= 0. - - // (x) : (exists y = [x / 2], z : x = 5y). - FlatAffineConstraints fac2(1); - fac2.addLocalFloorDiv({1, 0}, 2); // y = [x / 2]. - fac2.addEquality({1, -5, 0}); // x = 5y. - fac2.appendLocalId(); // Add local id z. - - fac1.mergeLocalIds(fac2); - - // Local space should be same. - EXPECT_EQ(fac1.getNumLocalIds(), fac2.getNumLocalIds()); - - // 1 division should be matched + 2 unmatched local ids. - EXPECT_EQ(fac1.getNumLocalIds(), 3u); - EXPECT_EQ(fac2.getNumLocalIds(), 3u); - } - - { - // (x) : (exists z = [x / 5], y = [x / 2] : x = 3y). - FlatAffineConstraints fac1(1); - fac1.addLocalFloorDiv({1, 0}, 5); // z = [x / 5]. - fac1.addLocalFloorDiv({1, 0, 0}, 2); // y = [x / 2]. - fac1.addEquality({1, 0, -3, 0}); // x = 3y. - - // (x) : (exists y = [x / 2], z = [x / 5]: x = 5z). - FlatAffineConstraints fac2(1); - fac2.addLocalFloorDiv({1, 0}, 2); // y = [x / 2]. - fac2.addLocalFloorDiv({1, 0, 0}, 5); // z = [x / 5]. - fac2.addEquality({1, 0, -5, 0}); // x = 5z. - - fac1.mergeLocalIds(fac2); - - // Local space should be same. - EXPECT_EQ(fac1.getNumLocalIds(), fac2.getNumLocalIds()); - - // 2 divisions should be matched. - EXPECT_EQ(fac1.getNumLocalIds(), 2u); - EXPECT_EQ(fac2.getNumLocalIds(), 2u); - } -} - -TEST(FlatAffineConstraintsTest, mergeDivisionsNestedDivsions) { - { - // (x) : (exists y = [x / 2], z = [x + y / 3]: y + z >= x). - FlatAffineConstraints fac1(1); - fac1.addLocalFloorDiv({1, 0}, 2); // y = [x / 2]. - fac1.addLocalFloorDiv({1, 1, 0}, 3); // z = [x + y / 3]. - fac1.addInequality({-1, 1, 1, 0}); // y + z >= x. - - // (x) : (exists y = [x / 2], z = [x + y / 3]: y + z <= x). - FlatAffineConstraints fac2(1); - fac2.addLocalFloorDiv({1, 0}, 2); // y = [x / 2]. - fac2.addLocalFloorDiv({1, 1, 0}, 3); // z = [x + y / 3]. - fac2.addInequality({1, -1, -1, 0}); // y + z <= x. - - fac1.mergeLocalIds(fac2); - - // Local space should be same. - EXPECT_EQ(fac1.getNumLocalIds(), fac2.getNumLocalIds()); - - // 2 divisions should be matched. - EXPECT_EQ(fac1.getNumLocalIds(), 2u); - EXPECT_EQ(fac2.getNumLocalIds(), 2u); - } - - { - // (x) : (exists y = [x / 2], z = [x + y / 3], w = [z + 1 / 5]: y + z >= x). - FlatAffineConstraints fac1(1); - fac1.addLocalFloorDiv({1, 0}, 2); // y = [x / 2]. - fac1.addLocalFloorDiv({1, 1, 0}, 3); // z = [x + y / 3]. - fac1.addLocalFloorDiv({0, 0, 1, 1}, 5); // w = [z + 1 / 5]. - fac1.addInequality({-1, 1, 1, 0, 0}); // y + z >= x. - - // (x) : (exists y = [x / 2], z = [x + y / 3], w = [z + 1 / 5]: y + z <= x). - FlatAffineConstraints fac2(1); - fac2.addLocalFloorDiv({1, 0}, 2); // y = [x / 2]. - fac2.addLocalFloorDiv({1, 1, 0}, 3); // z = [x + y / 3]. - fac2.addLocalFloorDiv({0, 0, 1, 1}, 5); // w = [z + 1 / 5]. - fac2.addInequality({1, -1, -1, 0, 0}); // y + z <= x. - - fac1.mergeLocalIds(fac2); - - // Local space should be same. - EXPECT_EQ(fac1.getNumLocalIds(), fac2.getNumLocalIds()); - - // 3 divisions should be matched. - EXPECT_EQ(fac1.getNumLocalIds(), 3u); - EXPECT_EQ(fac2.getNumLocalIds(), 3u); - } -} - -TEST(FlatAffineConstraintsTest, mergeDivisionsConstants) { - { - // (x) : (exists y = [x + 1 / 3], z = [x + 2 / 3]: y + z >= x). - FlatAffineConstraints fac1(1); - fac1.addLocalFloorDiv({1, 1}, 2); // y = [x + 1 / 2]. - fac1.addLocalFloorDiv({1, 0, 2}, 3); // z = [x + 2 / 3]. - fac1.addInequality({-1, 1, 1, 0}); // y + z >= x. - - // (x) : (exists y = [x + 1 / 3], z = [x + 2 / 3]: y + z <= x). - FlatAffineConstraints fac2(1); - fac2.addLocalFloorDiv({1, 1}, 2); // y = [x + 1 / 2]. - fac2.addLocalFloorDiv({1, 0, 2}, 3); // z = [x + 2 / 3]. - fac2.addInequality({1, -1, -1, 0}); // y + z <= x. - - fac1.mergeLocalIds(fac2); - - // Local space should be same. - EXPECT_EQ(fac1.getNumLocalIds(), fac2.getNumLocalIds()); - - // 2 divisions should be matched. - EXPECT_EQ(fac1.getNumLocalIds(), 2u); - EXPECT_EQ(fac2.getNumLocalIds(), 2u); - } -} - -} // namespace mlir diff --git a/mlir/unittests/Analysis/CMakeLists.txt b/mlir/unittests/Analysis/CMakeLists.txt --- a/mlir/unittests/Analysis/CMakeLists.txt +++ b/mlir/unittests/Analysis/CMakeLists.txt @@ -1,7 +1,6 @@ add_mlir_unittest(MLIRAnalysisTests AffineStructuresParser.cpp AffineStructuresParserTest.cpp - AffineStructuresTest.cpp PresburgerSetTest.cpp ) diff --git a/mlir/unittests/Analysis/Presburger/IntegerPolyhedronTest.cpp b/mlir/unittests/Analysis/Presburger/IntegerPolyhedronTest.cpp --- a/mlir/unittests/Analysis/Presburger/IntegerPolyhedronTest.cpp +++ b/mlir/unittests/Analysis/Presburger/IntegerPolyhedronTest.cpp @@ -7,14 +7,20 @@ //===----------------------------------------------------------------------===// #include "mlir/Analysis/Presburger/IntegerPolyhedron.h" +#include "../AffineStructuresParser.h" +#include "mlir/IR/MLIRContext.h" #include #include +#include + namespace mlir { using testing::ElementsAre; +enum class TestFunction { Sample, Empty }; + /// Construct a IntegerPolyhedron from a set of inequality and /// equality constraints. static IntegerPolyhedron @@ -30,6 +36,79 @@ return set; } +/// If fn is TestFunction::Sample (default): +/// If hasSample is true, check that findIntegerSample returns a valid sample +/// for the IntegerPolyhedron poly. +/// If hasSample is false, check that findIntegerSample returns None. +/// +/// If fn is TestFunction::Empty, check that isIntegerEmpty returns the +/// opposite of hasSample. +static void checkSample(bool hasSample, const IntegerPolyhedron &poly, + TestFunction fn = TestFunction::Sample) { + Optional> maybeSample; + switch (fn) { + case TestFunction::Sample: + maybeSample = poly.findIntegerSample(); + if (!hasSample) { + EXPECT_FALSE(maybeSample.hasValue()); + if (maybeSample.hasValue()) { + for (auto x : *maybeSample) + llvm::errs() << x << ' '; + llvm::errs() << '\n'; + } + } else { + ASSERT_TRUE(maybeSample.hasValue()); + EXPECT_TRUE(poly.containsPoint(*maybeSample)); + } + break; + case TestFunction::Empty: + EXPECT_EQ(!hasSample, poly.isIntegerEmpty()); + break; + } +} + +/// Check sampling for all the permutations of the dimensions for the given +/// constraint set. Since the GBR algorithm progresses dimension-wise, different +/// orderings may cause the algorithm to proceed differently. At least some of +///.these permutations should make it past the heuristics and test the +/// implementation of the GBR algorithm itself. +/// Use TestFunction fn to test. +static void checkPermutationsSample(bool hasSample, unsigned nDim, + ArrayRef> ineqs, + ArrayRef> eqs, + TestFunction fn = TestFunction::Sample) { + SmallVector perm(nDim); + std::iota(perm.begin(), perm.end(), 0); + auto permute = [&perm](ArrayRef coeffs) { + SmallVector permuted; + for (unsigned id : perm) + permuted.push_back(coeffs[id]); + permuted.push_back(coeffs.back()); + return permuted; + }; + do { + SmallVector, 4> permutedIneqs, permutedEqs; + for (const auto &ineq : ineqs) + permutedIneqs.push_back(permute(ineq)); + for (const auto &eq : eqs) + permutedEqs.push_back(permute(eq)); + + checkSample(hasSample, + makeSetFromConstraints(nDim, permutedIneqs, permutedEqs), fn); + } while (std::next_permutation(perm.begin(), perm.end())); +} + +/// Parses a IntegerPolyhedron from a StringRef. It is expected that the +/// string represents a valid IntegerSet, otherwise it will violate a gtest +/// assertion. +static IntegerPolyhedron parsePoly(StringRef str, MLIRContext *context) { + FailureOr poly = parseIntegerSetToFAC(str, context); + + EXPECT_TRUE(succeeded(poly)); + + return *poly; +} + TEST(IntegerPolyhedronTest, removeInequality) { IntegerPolyhedron set = makeSetFromConstraints(1, {{0, 0}, {1, 1}, {2, 2}, {3, 3}, {4, 4}}, {}); @@ -100,4 +179,722 @@ EXPECT_THAT(set.getInequality(0), testing::ElementsAre(12, 20, 40)); } +TEST(IntegerPolyhedronTest, FindSampleTest) { + // Bounded sets with only inequalities. + + MLIRContext context; + + // 0 <= 7x <= 5 + checkSample(true, parsePoly("(x) : (7 * x >= 0, -7 * x + 5 >= 0)", &context)); + + // 1 <= 5x and 5x <= 4 (no solution). + checkSample(false, + parsePoly("(x) : (5 * x - 1 >= 0, -5 * x + 4 >= 0)", &context)); + + // 1 <= 5x and 5x <= 9 (solution: x = 1). + checkSample(true, + parsePoly("(x) : (5 * x - 1 >= 0, -5 * x + 9 >= 0)", &context)); + + // Bounded sets with equalities. + // x >= 8 and 40 >= y and x = y. + checkSample(true, parsePoly("(x,y) : (x - 8 >= 0, -y + 40 >= 0, x - y == 0)", + &context)); + + // x <= 10 and y <= 10 and 10 <= z and x + 2y = 3z. + // solution: x = y = z = 10. + checkSample(true, parsePoly("(x,y,z) : (-x + 10 >= 0, -y + 10 >= 0, " + "z - 10 >= 0, x + 2 * y - 3 * z == 0)", + &context)); + + // x <= 10 and y <= 10 and 11 <= z and x + 2y = 3z. + // This implies x + 2y >= 33 and x + 2y <= 30, which has no solution. + checkSample(false, parsePoly("(x,y,z) : (-x + 10 >= 0, -y + 10 >= 0, " + "z - 11 >= 0, x + 2 * y - 3 * z == 0)", + &context)); + + // 0 <= r and r <= 3 and 4q + r = 7. + // Solution: q = 1, r = 3. + checkSample( + true, + parsePoly("(q,r) : (r >= 0, -r + 3 >= 0, 4 * q + r - 7 == 0)", &context)); + + // 4q + r = 7 and r = 0. + // Solution: q = 1, r = 3. + checkSample(false, + parsePoly("(q,r) : (4 * q + r - 7 == 0, r == 0)", &context)); + + // The next two sets are large sets that should take a long time to sample + // with a naive branch and bound algorithm but can be sampled efficiently with + // the GBR algorithm. + // + // This is a triangle with vertices at (1/3, 0), (2/3, 0) and (10000, 10000). + checkSample(true, parsePoly("(x,y) : (y >= 0, " + "300000 * x - 299999 * y - 100000 >= 0, " + "-300000 * x + 299998 * y + 200000 >= 0)", + &context)); + + // This is a tetrahedron with vertices at + // (1/3, 0, 0), (2/3, 0, 0), (2/3, 0, 10000), and (10000, 10000, 10000). + // The first three points form a triangular base on the xz plane with the + // apex at the fourth point, which is the only integer point. + checkPermutationsSample( + true, 3, + { + {0, 1, 0, 0}, // y >= 0 + {0, -1, 1, 0}, // z >= y + {300000, -299998, -1, + -100000}, // -300000x + 299998y + 100000 + z <= 0. + {-150000, 149999, 0, 100000}, // -150000x + 149999y + 100000 >= 0. + }, + {}); + + // Same thing with some spurious extra dimensions equated to constants. + checkSample( + true, + parsePoly("(a,b,c,d,e) : (b + d - e >= 0, -b + c - d + e >= 0, " + "300000 * a - 299998 * b - c - 9 * d + 21 * e - 112000 >= 0, " + "-150000 * a + 149999 * b - 15 * d + 47 * e + 68000 >= 0, " + "d - e == 0, d + e - 2000 == 0)", + &context)); + + // This is a tetrahedron with vertices at + // (1/3, 0, 0), (2/3, 0, 0), (2/3, 0, 100), (100, 100 - 1/3, 100). + checkPermutationsSample(false, 3, + { + {0, 1, 0, 0}, + {0, -300, 299, 0}, + {300 * 299, -89400, -299, -100 * 299}, + {-897, 894, 0, 598}, + }, + {}); + + // Two tests involving equalities that are integer empty but not rational + // empty. + + // This is a line segment from (0, 1/3) to (100, 100 + 1/3). + checkSample( + false, + parsePoly("(x,y) : (x >= 0, -x + 100 >= 0, 3 * x - 3 * y + 1 == 0)", + &context)); + + // A thin parallelogram. 0 <= x <= 100 and x + 1/3 <= y <= x + 2/3. + checkSample(false, + parsePoly("(x,y) : (x >= 0, -x + 100 >= 0, " + "3 * x - 3 * y + 2 >= 0, -3 * x + 3 * y - 1 >= 0)", + &context)); + + checkSample(true, parsePoly("(x,y) : (2 * x >= 0, -2 * x + 99 >= 0, " + "2 * y >= 0, -2 * y + 99 >= 0)", + &context)); + + // 2D cone with apex at (10000, 10000) and + // edges passing through (1/3, 0) and (2/3, 0). + checkSample(true, parsePoly("(x,y) : (300000 * x - 299999 * y - 100000 >= 0, " + "-300000 * x + 299998 * y + 200000 >= 0)", + &context)); + + // Cartesian product of a tetrahedron and a 2D cone. + // The tetrahedron has vertices at + // (1/3, 0, 0), (2/3, 0, 0), (2/3, 0, 10000), and (10000, 10000, 10000). + // The first three points form a triangular base on the xz plane with the + // apex at the fourth point, which is the only integer point. + // The cone has apex at (10000, 10000) and + // edges passing through (1/3, 0) and (2/3, 0). + checkPermutationsSample( + true /* not empty */, 5, + { + // Tetrahedron contraints: + {0, 1, 0, 0, 0, 0}, // y >= 0 + {0, -1, 1, 0, 0, 0}, // z >= y + // -300000x + 299998y + 100000 + z <= 0. + {300000, -299998, -1, 0, 0, -100000}, + // -150000x + 149999y + 100000 >= 0. + {-150000, 149999, 0, 0, 0, 100000}, + + // Triangle constraints: + // 300000p - 299999q >= 100000 + {0, 0, 0, 300000, -299999, -100000}, + // -300000p + 299998q + 200000 >= 0 + {0, 0, 0, -300000, 299998, 200000}, + }, + {}); + + // Cartesian product of same tetrahedron as above and {(p, q) : 1/3 <= p <= + // 2/3}. Since the second set is empty, the whole set is too. + checkPermutationsSample( + false /* empty */, 5, + { + // Tetrahedron contraints: + {0, 1, 0, 0, 0, 0}, // y >= 0 + {0, -1, 1, 0, 0, 0}, // z >= y + // -300000x + 299998y + 100000 + z <= 0. + {300000, -299998, -1, 0, 0, -100000}, + // -150000x + 149999y + 100000 >= 0. + {-150000, 149999, 0, 0, 0, 100000}, + + // Second set constraints: + // 3p >= 1 + {0, 0, 0, 3, 0, -1}, + // 3p <= 2 + {0, 0, 0, -3, 0, 2}, + }, + {}); + + // Cartesian product of same tetrahedron as above and + // {(p, q, r) : 1 <= p <= 2 and p = 3q + 3r}. + // Since the second set is empty, the whole set is too. + checkPermutationsSample( + false /* empty */, 5, + { + // Tetrahedron contraints: + {0, 1, 0, 0, 0, 0, 0}, // y >= 0 + {0, -1, 1, 0, 0, 0, 0}, // z >= y + // -300000x + 299998y + 100000 + z <= 0. + {300000, -299998, -1, 0, 0, 0, -100000}, + // -150000x + 149999y + 100000 >= 0. + {-150000, 149999, 0, 0, 0, 0, 100000}, + + // Second set constraints: + // p >= 1 + {0, 0, 0, 1, 0, 0, -1}, + // p <= 2 + {0, 0, 0, -1, 0, 0, 2}, + }, + { + {0, 0, 0, 1, -3, -3, 0}, // p = 3q + 3r + }); + + // Cartesian product of a tetrahedron and a 2D cone. + // The tetrahedron is empty and has vertices at + // (1/3, 0, 0), (2/3, 0, 0), (2/3, 0, 100), and (100, 100 - 1/3, 100). + // The cone has apex at (10000, 10000) and + // edges passing through (1/3, 0) and (2/3, 0). + // Since the tetrahedron is empty, the Cartesian product is too. + checkPermutationsSample(false /* empty */, 5, + { + // Tetrahedron contraints: + {0, 1, 0, 0, 0, 0}, + {0, -300, 299, 0, 0, 0}, + {300 * 299, -89400, -299, 0, 0, -100 * 299}, + {-897, 894, 0, 0, 0, 598}, + + // Triangle constraints: + // 300000p - 299999q >= 100000 + {0, 0, 0, 300000, -299999, -100000}, + // -300000p + 299998q + 200000 >= 0 + {0, 0, 0, -300000, 299998, 200000}, + }, + {}); + + // Cartesian product of same tetrahedron as above and + // {(p, q) : 1/3 <= p <= 2/3}. + checkPermutationsSample(false /* empty */, 5, + { + // Tetrahedron contraints: + {0, 1, 0, 0, 0, 0}, + {0, -300, 299, 0, 0, 0}, + {300 * 299, -89400, -299, 0, 0, -100 * 299}, + {-897, 894, 0, 0, 0, 598}, + + // Second set constraints: + // 3p >= 1 + {0, 0, 0, 3, 0, -1}, + // 3p <= 2 + {0, 0, 0, -3, 0, 2}, + }, + {}); + + checkSample(true, parsePoly("(x, y, z) : (2 * x - 1 >= 0, x - y - 1 == 0, " + "y - z == 0)", + &context)); + + // Regression tests for the computation of dual coefficients. + checkSample(false, parsePoly("(x, y, z) : (" + "6*x - 4*y + 9*z + 2 >= 0," + "x + 5*y + z + 5 >= 0," + "-4*x + y + 2*z - 1 >= 0," + "-3*x - 2*y - 7*z - 1 >= 0," + "-7*x - 5*y - 9*z - 1 >= 0)", + &context)); + checkSample(true, parsePoly("(x, y, z) : (" + "3*x + 3*y + 3 >= 0," + "-4*x - 8*y - z + 4 >= 0," + "-7*x - 4*y + z + 1 >= 0," + "2*x - 7*y - 8*z - 7 >= 0," + "9*x + 8*y - 9*z - 7 >= 0)", + &context)); +} + +TEST(IntegerPolyhedronTest, IsIntegerEmptyTest) { + + MLIRContext context; + + // 1 <= 5x and 5x <= 4 (no solution). + EXPECT_TRUE(parsePoly("(x) : (5 * x - 1 >= 0, -5 * x + 4 >= 0)", &context) + .isIntegerEmpty()); + // 1 <= 5x and 5x <= 9 (solution: x = 1). + EXPECT_FALSE(parsePoly("(x) : (5 * x - 1 >= 0, -5 * x + 9 >= 0)", &context) + .isIntegerEmpty()); + + // Unbounded sets. + EXPECT_TRUE(parsePoly("(x,y,z) : (2 * y - 1 >= 0, -2 * y + 1 >= 0, " + "2 * z - 1 >= 0, 2 * x - 1 == 0)", + &context) + .isIntegerEmpty()); + + EXPECT_FALSE(parsePoly("(x,y,z) : (2 * x - 1 >= 0, -3 * x + 3 >= 0, " + "5 * z - 6 >= 0, -7 * z + 17 >= 0, 3 * y - 2 >= 0)", + &context) + .isIntegerEmpty()); + + EXPECT_FALSE( + parsePoly("(x,y,z) : (2 * x - 1 >= 0, x - y - 1 == 0, y - z == 0)", + &context) + .isIntegerEmpty()); + + // IntegerPolyhedron::isEmpty() does not detect the following sets to be + // empty. + + // 3x + 7y = 1 and 0 <= x, y <= 10. + // Since x and y are non-negative, 3x + 7y can never be 1. + EXPECT_TRUE(parsePoly("(x,y) : (x >= 0, -x + 10 >= 0, y >= 0, -y + 10 >= 0, " + "3 * x + 7 * y - 1 == 0)", + &context) + .isIntegerEmpty()); + + // 2x = 3y and y = x - 1 and x + y = 6z + 2 and 0 <= x, y <= 100. + // Substituting y = x - 1 in 3y = 2x, we obtain x = 3 and hence y = 2. + // Since x + y = 5 cannot be equal to 6z + 2 for any z, the set is empty. + EXPECT_TRUE( + parsePoly("(x,y,z) : (x >= 0, -x + 100 >= 0, y >= 0, -y + 100 >= 0, " + "2 * x - 3 * y == 0, x - y - 1 == 0, x + y - 6 * z - 2 == 0)", + &context) + .isIntegerEmpty()); + + // 2x = 3y and y = x - 1 + 6z and x + y = 6q + 2 and 0 <= x, y <= 100. + // 2x = 3y implies x is a multiple of 3 and y is even. + // Now y = x - 1 + 6z implies y = 2 mod 3. In fact, since y is even, we have + // y = 2 mod 6. Then since x = y + 1 + 6z, we have x = 3 mod 6, implying + // x + y = 5 mod 6, which contradicts x + y = 6q + 2, so the set is empty. + EXPECT_TRUE( + parsePoly( + "(x,y,z,q) : (x >= 0, -x + 100 >= 0, y >= 0, -y + 100 >= 0, " + "2 * x - 3 * y == 0, x - y + 6 * z - 1 == 0, x + y - 6 * q - 2 == 0)", + &context) + .isIntegerEmpty()); + + // Set with symbols. + EXPECT_FALSE(parsePoly("(x)[s] : (x + s >= 0, x - s == 0)", &context) + .isIntegerEmpty()); +} + +TEST(IntegerPolyhedronTest, removeRedundantConstraintsTest) { + MLIRContext context; + + IntegerPolyhedron poly = + parsePoly("(x) : (x - 2 >= 0, -x + 2 >= 0, x - 2 == 0)", &context); + poly.removeRedundantConstraints(); + + // Both inequalities are redundant given the equality. Both have been removed. + EXPECT_EQ(poly.getNumInequalities(), 0u); + EXPECT_EQ(poly.getNumEqualities(), 1u); + + IntegerPolyhedron poly2 = + parsePoly("(x,y) : (x - 3 >= 0, y - 2 >= 0, x - y == 0)", &context); + poly2.removeRedundantConstraints(); + + // The second inequality is redundant and should have been removed. The + // remaining inequality should be the first one. + EXPECT_EQ(poly2.getNumInequalities(), 1u); + EXPECT_THAT(poly2.getInequality(0), ElementsAre(1, 0, -3)); + EXPECT_EQ(poly2.getNumEqualities(), 1u); + + IntegerPolyhedron poly3 = + parsePoly("(x,y,z) : (x - y == 0, x - z == 0, y - z == 0)", &context); + poly3.removeRedundantConstraints(); + + // One of the three equalities can be removed. + EXPECT_EQ(poly3.getNumInequalities(), 0u); + EXPECT_EQ(poly3.getNumEqualities(), 2u); + + IntegerPolyhedron poly4 = + parsePoly("(a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q) : (" + "b - 1 >= 0," + "-b + 500 >= 0," + "-16 * d + f >= 0," + "f - 1 >= 0," + "-f + 998 >= 0," + "16 * d - f + 15 >= 0," + "-16 * e + g >= 0," + "g - 1 >= 0," + "-g + 998 >= 0," + "16 * e - g + 15 >= 0," + "h >= 0," + "-h + 1 >= 0," + "j - 1 >= 0," + "-j + 500 >= 0," + "-f + 16 * l + 15 >= 0," + "f - 16 * l >= 0," + "-16 * m + o >= 0," + "o - 1 >= 0," + "-o + 998 >= 0," + "16 * m - o + 15 >= 0," + "p >= 0," + "-p + 1 >= 0," + "-g - h + 8 * q + 8 >= 0," + "-o - p + 8 * q + 8 >= 0," + "o + p - 8 * q - 1 >= 0," + "g + h - 8 * q - 1 >= 0," + "-f + n >= 0," + "f - n >= 0," + "k - 10 >= 0," + "-k + 10 >= 0," + "i - 13 >= 0," + "-i + 13 >= 0," + "c - 10 >= 0," + "-c + 10 >= 0," + "a - 13 >= 0," + "-a + 13 >= 0" + ")", + &context); + + // The above is a large set of constraints without any redundant constraints, + // as verified by the Fourier-Motzkin based removeRedundantInequalities. + unsigned nIneq = poly4.getNumInequalities(); + unsigned nEq = poly4.getNumEqualities(); + poly4.removeRedundantInequalities(); + ASSERT_EQ(poly4.getNumInequalities(), nIneq); + ASSERT_EQ(poly4.getNumEqualities(), nEq); + // Now we test that removeRedundantConstraints does not find any constraints + // to be redundant either. + poly4.removeRedundantConstraints(); + EXPECT_EQ(poly4.getNumInequalities(), nIneq); + EXPECT_EQ(poly4.getNumEqualities(), nEq); + + IntegerPolyhedron poly5 = parsePoly( + "(x,y) : (128 * x + 127 >= 0, -x + 7 >= 0, -128 * x + y >= 0, y >= 0)", + &context); + // 128x + 127 >= 0 implies that 128x >= 0, since x has to be an integer. + // (This should be caught by GCDTightenInqualities().) + // So -128x + y >= 0 and 128x + 127 >= 0 imply y >= 0 since we have + // y >= 128x >= 0. + poly5.removeRedundantConstraints(); + EXPECT_EQ(poly5.getNumInequalities(), 3u); + SmallVector redundantConstraint = {0, 1, 0}; + for (unsigned i = 0; i < 3; ++i) { + // Ensure that the removed constraint was the redundant constraint [3]. + EXPECT_NE(poly5.getInequality(i), ArrayRef(redundantConstraint)); + } +} + +TEST(IntegerPolyhedronTest, addConstantUpperBound) { + IntegerPolyhedron poly(2); + poly.addBound(IntegerPolyhedron::UB, 0, 1); + EXPECT_EQ(poly.atIneq(0, 0), -1); + EXPECT_EQ(poly.atIneq(0, 1), 0); + EXPECT_EQ(poly.atIneq(0, 2), 1); + + poly.addBound(IntegerPolyhedron::UB, {1, 2, 3}, 1); + EXPECT_EQ(poly.atIneq(1, 0), -1); + EXPECT_EQ(poly.atIneq(1, 1), -2); + EXPECT_EQ(poly.atIneq(1, 2), -2); +} + +TEST(IntegerPolyhedronTest, addConstantLowerBound) { + IntegerPolyhedron poly(2); + poly.addBound(IntegerPolyhedron::LB, 0, 1); + EXPECT_EQ(poly.atIneq(0, 0), 1); + EXPECT_EQ(poly.atIneq(0, 1), 0); + EXPECT_EQ(poly.atIneq(0, 2), -1); + + poly.addBound(IntegerPolyhedron::LB, {1, 2, 3}, 1); + EXPECT_EQ(poly.atIneq(1, 0), 1); + EXPECT_EQ(poly.atIneq(1, 1), 2); + EXPECT_EQ(poly.atIneq(1, 2), 2); +} + +/// Check if the expected division representation of local variables matches the +/// computed representation. The expected division representation is given as +/// a vector of expressions set in `expectedDividends` and the corressponding +/// denominator in `expectedDenominators`. The `denominators` and `dividends` +/// obtained through `getLocalRepr` function is verified against the +/// `expectedDenominators` and `expectedDividends` respectively. +static void checkDivisionRepresentation( + IntegerPolyhedron &poly, + const std::vector> &expectedDividends, + const SmallVectorImpl &expectedDenominators) { + + std::vector> dividends; + SmallVector denominators; + + poly.getLocalReprs(dividends, denominators); + + // Check that the `denominators` and `expectedDenominators` match. + EXPECT_TRUE(expectedDenominators == denominators); + + // Check that the `dividends` and `expectedDividends` match. If the + // denominator for a division is zero, we ignore its dividend. + EXPECT_TRUE(dividends.size() == expectedDividends.size()); + for (unsigned i = 0, e = dividends.size(); i < e; ++i) + if (denominators[i] != 0) + EXPECT_TRUE(expectedDividends[i] == dividends[i]); +} + +TEST(IntegerPolyhedronTest, computeLocalReprSimple) { + IntegerPolyhedron poly(1); + + poly.addLocalFloorDiv({1, 4}, 10); + poly.addLocalFloorDiv({1, 0, 100}, 10); + + std::vector> divisions = {{1, 0, 0, 4}, + {1, 0, 0, 100}}; + SmallVector denoms = {10, 10}; + + // Check if floordivs can be computed when no other inequalities exist + // and floor divs do not depend on each other. + checkDivisionRepresentation(poly, divisions, denoms); +} + +TEST(IntegerPolyhedronTest, computeLocalReprConstantFloorDiv) { + IntegerPolyhedron poly(4); + + poly.addInequality({1, 0, 3, 1, 2}); + poly.addInequality({1, 2, -8, 1, 10}); + poly.addEquality({1, 2, -4, 1, 10}); + + poly.addLocalFloorDiv({0, 0, 0, 0, 10}, 30); + poly.addLocalFloorDiv({0, 0, 0, 0, 0, 99}, 101); + + std::vector> divisions = {{0, 0, 0, 0, 0, 0, 10}, + {0, 0, 0, 0, 0, 0, 99}}; + SmallVector denoms = {30, 101}; + + // Check if floordivs with constant numerator can be computed. + checkDivisionRepresentation(poly, divisions, denoms); +} + +TEST(IntegerPolyhedronTest, computeLocalReprRecursive) { + IntegerPolyhedron poly(4); + poly.addInequality({1, 0, 3, 1, 2}); + poly.addInequality({1, 2, -8, 1, 10}); + poly.addEquality({1, 2, -4, 1, 10}); + + poly.addLocalFloorDiv({0, -2, 7, 2, 10}, 3); + poly.addLocalFloorDiv({3, 0, 9, 2, 2, 10}, 5); + poly.addLocalFloorDiv({0, 1, -123, 2, 0, -4, 10}, 3); + + poly.addInequality({1, 2, -2, 1, -5, 0, 6, 100}); + poly.addInequality({1, 2, -8, 1, 3, 7, 0, -9}); + + std::vector> divisions = { + {0, -2, 7, 2, 0, 0, 0, 10}, + {3, 0, 9, 2, 2, 0, 0, 10}, + {0, 1, -123, 2, 0, -4, 0, 10}}; + + SmallVector denoms = {3, 5, 3}; + + // Check if floordivs which may depend on other floordivs can be computed. + checkDivisionRepresentation(poly, divisions, denoms); +} + +TEST(IntegerPolyhedronTest, computeLocalReprTightUpperBound) { + MLIRContext context; + + { + IntegerPolyhedron poly = parsePoly("(i) : (i mod 3 - 1 >= 0)", &context); + + // The set formed by the poly is: + // 3q - i + 2 >= 0 <-- Division lower bound + // -3q + i - 1 >= 0 + // -3q + i >= 0 <-- Division upper bound + // We remove redundant constraints to get the set: + // 3q - i + 2 >= 0 <-- Division lower bound + // -3q + i - 1 >= 0 <-- Tighter division upper bound + // thus, making the upper bound tighter. + poly.removeRedundantConstraints(); + + std::vector> divisions = {{1, 0, 0}}; + SmallVector denoms = {3}; + + // Check if the divisions can be computed even with a tighter upper bound. + checkDivisionRepresentation(poly, divisions, denoms); + } + + { + IntegerPolyhedron poly = parsePoly( + "(i, j, q) : (4*q - i - j + 2 >= 0, -4*q + i + j >= 0)", &context); + // Convert `q` to a local variable. + poly.convertDimToLocal(2, 3); + + std::vector> divisions = {{1, 1, 0, 1}}; + SmallVector denoms = {4}; + + // Check if the divisions can be computed even with a tighter upper bound. + checkDivisionRepresentation(poly, divisions, denoms); + } +} + +TEST(IntegerPolyhedronTest, computeLocalReprNoRepr) { + MLIRContext context; + IntegerPolyhedron poly = + parsePoly("(x, q) : (x - 3 * q >= 0, -x + 3 * q + 3 >= 0)", &context); + // Convert q to a local variable. + poly.convertDimToLocal(1, 2); + + std::vector> divisions = {{0, 0, 0}}; + SmallVector denoms = {0}; + + // Check that no division is computed. + checkDivisionRepresentation(poly, divisions, denoms); +} + +TEST(IntegerPolyhedronTest, simplifyLocalsTest) { + // (x) : (exists y: 2x + y = 1 and y = 2). + IntegerPolyhedron poly(1, 0, 1); + poly.addEquality({2, 1, -1}); + poly.addEquality({0, 1, -2}); + + EXPECT_TRUE(poly.isEmpty()); + + // (x) : (exists y, z, w: 3x + y = 1 and 2y = z and 3y = w and z = w). + IntegerPolyhedron poly2(1, 0, 3); + poly2.addEquality({3, 1, 0, 0, -1}); + poly2.addEquality({0, 2, -1, 0, 0}); + poly2.addEquality({0, 3, 0, -1, 0}); + poly2.addEquality({0, 0, 1, -1, 0}); + + EXPECT_TRUE(poly2.isEmpty()); + + // (x) : (exists y: x >= y + 1 and 2x + y = 0 and y >= -1). + IntegerPolyhedron poly3(1, 0, 1); + poly3.addInequality({1, -1, -1}); + poly3.addInequality({0, 1, 1}); + poly3.addEquality({2, 1, 0}); + + EXPECT_TRUE(poly3.isEmpty()); +} + +TEST(IntegerPolyhedronTest, mergeDivisionsSimple) { + { + // (x) : (exists z, y = [x / 2] : x = 3y and x + z + 1 >= 0). + IntegerPolyhedron poly1(1, 0, 1); + poly1.addLocalFloorDiv({1, 0, 0}, 2); // y = [x / 2]. + poly1.addEquality({1, 0, -3, 0}); // x = 3y. + poly1.addInequality({1, 1, 0, 1}); // x + z + 1 >= 0. + + // (x) : (exists y = [x / 2], z : x = 5y). + IntegerPolyhedron poly2(1); + poly2.addLocalFloorDiv({1, 0}, 2); // y = [x / 2]. + poly2.addEquality({1, -5, 0}); // x = 5y. + poly2.appendLocalId(); // Add local id z. + + poly1.mergeLocalIds(poly2); + + // Local space should be same. + EXPECT_EQ(poly1.getNumLocalIds(), poly2.getNumLocalIds()); + + // 1 division should be matched + 2 unmatched local ids. + EXPECT_EQ(poly1.getNumLocalIds(), 3u); + EXPECT_EQ(poly2.getNumLocalIds(), 3u); + } + + { + // (x) : (exists z = [x / 5], y = [x / 2] : x = 3y). + IntegerPolyhedron poly1(1); + poly1.addLocalFloorDiv({1, 0}, 5); // z = [x / 5]. + poly1.addLocalFloorDiv({1, 0, 0}, 2); // y = [x / 2]. + poly1.addEquality({1, 0, -3, 0}); // x = 3y. + + // (x) : (exists y = [x / 2], z = [x / 5]: x = 5z). + IntegerPolyhedron poly2(1); + poly2.addLocalFloorDiv({1, 0}, 2); // y = [x / 2]. + poly2.addLocalFloorDiv({1, 0, 0}, 5); // z = [x / 5]. + poly2.addEquality({1, 0, -5, 0}); // x = 5z. + + poly1.mergeLocalIds(poly2); + + // Local space should be same. + EXPECT_EQ(poly1.getNumLocalIds(), poly2.getNumLocalIds()); + + // 2 divisions should be matched. + EXPECT_EQ(poly1.getNumLocalIds(), 2u); + EXPECT_EQ(poly2.getNumLocalIds(), 2u); + } +} + +TEST(IntegerPolyhedronTest, mergeDivisionsNestedDivsions) { + { + // (x) : (exists y = [x / 2], z = [x + y / 3]: y + z >= x). + IntegerPolyhedron poly1(1); + poly1.addLocalFloorDiv({1, 0}, 2); // y = [x / 2]. + poly1.addLocalFloorDiv({1, 1, 0}, 3); // z = [x + y / 3]. + poly1.addInequality({-1, 1, 1, 0}); // y + z >= x. + + // (x) : (exists y = [x / 2], z = [x + y / 3]: y + z <= x). + IntegerPolyhedron poly2(1); + poly2.addLocalFloorDiv({1, 0}, 2); // y = [x / 2]. + poly2.addLocalFloorDiv({1, 1, 0}, 3); // z = [x + y / 3]. + poly2.addInequality({1, -1, -1, 0}); // y + z <= x. + + poly1.mergeLocalIds(poly2); + + // Local space should be same. + EXPECT_EQ(poly1.getNumLocalIds(), poly2.getNumLocalIds()); + + // 2 divisions should be matched. + EXPECT_EQ(poly1.getNumLocalIds(), 2u); + EXPECT_EQ(poly2.getNumLocalIds(), 2u); + } + + { + // (x) : (exists y = [x / 2], z = [x + y / 3], w = [z + 1 / 5]: y + z >= x). + IntegerPolyhedron poly1(1); + poly1.addLocalFloorDiv({1, 0}, 2); // y = [x / 2]. + poly1.addLocalFloorDiv({1, 1, 0}, 3); // z = [x + y / 3]. + poly1.addLocalFloorDiv({0, 0, 1, 1}, 5); // w = [z + 1 / 5]. + poly1.addInequality({-1, 1, 1, 0, 0}); // y + z >= x. + + // (x) : (exists y = [x / 2], z = [x + y / 3], w = [z + 1 / 5]: y + z <= x). + IntegerPolyhedron poly2(1); + poly2.addLocalFloorDiv({1, 0}, 2); // y = [x / 2]. + poly2.addLocalFloorDiv({1, 1, 0}, 3); // z = [x + y / 3]. + poly2.addLocalFloorDiv({0, 0, 1, 1}, 5); // w = [z + 1 / 5]. + poly2.addInequality({1, -1, -1, 0, 0}); // y + z <= x. + + poly1.mergeLocalIds(poly2); + + // Local space should be same. + EXPECT_EQ(poly1.getNumLocalIds(), poly2.getNumLocalIds()); + + // 3 divisions should be matched. + EXPECT_EQ(poly1.getNumLocalIds(), 3u); + EXPECT_EQ(poly2.getNumLocalIds(), 3u); + } +} + +TEST(IntegerPolyhedronTest, mergeDivisionsConstants) { + { + // (x) : (exists y = [x + 1 / 3], z = [x + 2 / 3]: y + z >= x). + IntegerPolyhedron poly1(1); + poly1.addLocalFloorDiv({1, 1}, 2); // y = [x + 1 / 2]. + poly1.addLocalFloorDiv({1, 0, 2}, 3); // z = [x + 2 / 3]. + poly1.addInequality({-1, 1, 1, 0}); // y + z >= x. + + // (x) : (exists y = [x + 1 / 3], z = [x + 2 / 3]: y + z <= x). + IntegerPolyhedron poly2(1); + poly2.addLocalFloorDiv({1, 1}, 2); // y = [x + 1 / 2]. + poly2.addLocalFloorDiv({1, 0, 2}, 3); // z = [x + 2 / 3]. + poly2.addInequality({1, -1, -1, 0}); // y + z <= x. + + poly1.mergeLocalIds(poly2); + + // Local space should be same. + EXPECT_EQ(poly1.getNumLocalIds(), poly2.getNumLocalIds()); + + // 2 divisions should be matched. + EXPECT_EQ(poly1.getNumLocalIds(), 2u); + EXPECT_EQ(poly2.getNumLocalIds(), 2u); + } +} + } // namespace mlir