Support emptiness checks for unbounded FlatAffineConstraints.

[arjunp]

Summary

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

[arjunp]

Added attribution for the algorithm.

[ftynse]

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

[arjunp]

Trying to fix the patch.

[arjunp]

Seems to be fixed now.

[arjunp]

Resolved clang-tidy issue.

[ftynse]

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

31 ↗

(On Diff #315321)

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

mlir/lib/Analysis/AffineStructures.cpp

1079

Shouldn't this be atEq?

1103

Let's not commit commented-out code

1106–1107

Or you can do

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

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

1185

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

mlir/lib/Analysis/Presburger/LinearTransform.cpp

14 ↗

(On Diff #315321)

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

16 ↗

(On Diff #315321)

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

23 ↗

(On Diff #315321)

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

23–24 ↗

(On Diff #315321)

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

75 ↗

(On Diff #315321)

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

110 ↗

(On Diff #315321)

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

mlir/unittests/Analysis/AffineStructuresTest.cpp

254 ↗

(On Diff #315321)

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

[arjunp]

Address comments.

[arjunp]

Thanks for the review!

mlir/lib/Analysis/AffineStructures.cpp

1079

Sorry, yes. Fixed and added a test case.

mlir/lib/Analysis/Presburger/LinearTransform.cpp

16 ↗

(On Diff #315321)

Yes, fixed.

23–24 ↗

(On Diff #315321)

Yeah, my intention was to use Euclidean division.

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

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

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

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

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

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

mlir/unittests/Analysis/AffineStructuresTest.cpp

254 ↗

(On Diff #315321)

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

[andydavis1]

mlir/lib/Analysis/AffineStructures.cpp

1065

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

[arjunp]

Addressed Andy's comment.

[bollu]

mlir/lib/Analysis/AffineStructures.cpp

1066

consider

const unsigned dirsNumCols = getNumCols() - 1;

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

mlir/lib/Analysis/Presburger/Matrix.cpp

88 ↗

(On Diff #316329)

nit: eliminate unnecessary return?

mlir/lib/Analysis/Presburger/Simplex.cpp

523 ↗

(On Diff #316329)

Consider:

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

instead of:

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

[arjunp]

Addressed Siddharth's commments.

[ftynse]

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

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

mlir/lib/Analysis/AffineStructures.cpp

1066

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

[arjunp]

Added unit tests for makeTransformToColumnEchelon.

[ftynse]

Thank you!

[arjunp]

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

mlir/lib/Analysis/AffineStructures.cpp

Matrix FlatAffineConstraints::getBoundedDirections() const {
   // The constructor adds the inequalities to the simplex first, so this
   // processes all the inequalities.
   for (unsigned i = 0, e = getNumInequalities(); i < e; ++i) {
     if (simplex.isBoundedAlongConstraint(i))
       boundedIneqs.push_back(i);
   }
 
   // The direction vector is given by the coefficients and does not include the
   // constant term, so the matrix has one fewer column.
   Matrix dirs(boundedIneqs.size() + getNumEqualities(), getNumCols() - 1);
 
   // Copy the bounded inequalities.
   unsigned row = 0;
   for (unsigned i : boundedIneqs) {
     for (unsigned col = 0, e = getNumCols() - 1; col < e; ++col)
       dirs(row, col) = atIneq(i, col);
     ++row;
   }
 
   // Copy the equalities. All the equalities' perpendiculars are bounded.
   for (unsigned i = 0, e = getNumEqualities(); i < e; ++i) {
     for (unsigned col = 0, e = getNumCols() - 1; col < e; ++col)
       dirs(row, col) = atIneq(i, col);
     ++row;
   }
 
   return dirs;
 }
 
 bool eqInvolvesSuffixDims(const FlatAffineConstraints &fac, unsigned eqIndex,
                           unsigned numDims) {
   for (unsigned e = fac.getNumDimIds(), j = e - numDims; j < e; ++j)
     if (fac.atEq(eqIndex, j) != 0)
       return true;
   return false;
 }
 bool ineqInvolvesSuffixDims(const FlatAffineConstraints &fac,
                             unsigned ineqIndex, unsigned numDims) {
   for (unsigned e = fac.getNumDimIds(), j = e - numDims; j < e; ++j)
     if (fac.atIneq(ineqIndex, j) != 0)
       return true;
   return false;
 }
 
 void removeConstraintsInvolvingSuffixDims(FlatAffineConstraints &fac,
                                           unsigned unboundedDims) {
   // unsigned unboundedDims = fac.getNumIds() - boundedDims;
 
   // We iterate backwards so that whether we remove constraint i or not, the
   // next constraint to be tested is always i - 1. i has to be signed to avoid
   // overflow at zero.
   for (int i = fac.getNumEqualities() - 1; i >= 0; i--)
     if (eqInvolvesSuffixDims(fac, i, unboundedDims))
       fac.removeEquality(i);
   for (int i = fac.getNumInequalities() - 1; i >= 0; i--)
     if (ineqInvolvesSuffixDims(fac, i, unboundedDims))
       fac.removeInequality(i);
 }
 
 /// Let this set be S. If S is bounded then we directly call into the GBR
-/// sampling algorithm. Otherwise, if there are some unbounded directions,
-/// vectors v such that S extends to infininty along v or -v, then a different
-/// algorithm is used:
+/// sampling algorithm. Otherwise, there are some unbounded directions, i.e.,
+/// vectors v such that S extends to infininty along v or -v. In this case we
+/// use an algorithm described in the integer set library (isl) manual and used
+/// by the isl_set_sample function in that library. The algorithm is:
 ///
 /// 1) Apply an invertible transform T to S to obtain S*T, such that all
 /// dimensions in which S*T is bounded lie in the linear span of a prefix of the
 /// dimensions.
 ///
 /// 2) Construct a set transformedSet by removing all constraints that involve
 /// the unbounded dimensions and also deleting the unbounded dimensions. Note
 /// that this is a bounded set.
bool FlatAffineConstraints::isIntegerEmpty() const {
   FlatAffineConstraints transformedSet = result.second.apply(*this);
 
   unsigned numBoundedDims = result.first;
   unsigned numUnboundedDims = getNumIds() - numBoundedDims;
   removeConstraintsInvolvingSuffixDims(transformedSet, numUnboundedDims);
   // Remove all the unbounded dimensions. i has to be signed to avoid overflow
   // if boundedDims is zero.
   for (int i = transformedSet.getNumIds() - 1; i >= int(numBoundedDims); i--)
     transformedSet.removeId(numBoundedDims);
 
   return !Simplex(transformedSet).findIntegerSample().hasValue();
 }
 
 Optional<SmallVector<int64_t, 8>>
 FlatAffineConstraints::findIntegerSample() const {
   return Simplex(*this).findIntegerSample();
 }