[Mlir-commits] [mlir] [MLIR][Presburger][WIP] Implement vertex enumeration and chamber decomposition for polytope generating function computation. (PR #78987)

llvmlistbot at llvm.org llvmlistbot at llvm.org
Wed Jan 24 09:31:36 PST 2024


https://github.com/Abhinav271828 updated https://github.com/llvm/llvm-project/pull/78987

>From be5fa415840591c7e001c18cbe0c3f2e0c9aa274 Mon Sep 17 00:00:00 2001
From: Abhinav271828 <abhinav.m at research.iiit.ac.in>
Date: Mon, 22 Jan 2024 14:56:23 +0530
Subject: [PATCH 01/26] Initial commit

---
 .../mlir/Analysis/Presburger/Barvinok.h       |   5 +
 .../Analysis/Presburger/GeneratingFunction.h  |   2 +-
 .../Analysis/Presburger/IntegerRelation.h     |   4 +
 mlir/lib/Analysis/Presburger/Barvinok.cpp     | 217 ++++++++++++++++++
 .../Analysis/Presburger/IntegerRelation.cpp   |  36 +++
 5 files changed, 263 insertions(+), 1 deletion(-)

diff --git a/mlir/include/mlir/Analysis/Presburger/Barvinok.h b/mlir/include/mlir/Analysis/Presburger/Barvinok.h
index b70ec33b693235f..6feb7e8f09233f6 100644
--- a/mlir/include/mlir/Analysis/Presburger/Barvinok.h
+++ b/mlir/include/mlir/Analysis/Presburger/Barvinok.h
@@ -27,6 +27,7 @@
 #include "mlir/Analysis/Presburger/GeneratingFunction.h"
 #include "mlir/Analysis/Presburger/IntegerRelation.h"
 #include "mlir/Analysis/Presburger/Matrix.h"
+#include "mlir/Analysis/Presburger/PresburgerRelation.h"
 #include "mlir/Analysis/Presburger/QuasiPolynomial.h"
 #include <optional>
 
@@ -84,6 +85,10 @@ ConeH getDual(ConeV cone);
 GeneratingFunction unimodularConeGeneratingFunction(ParamPoint vertex, int sign,
                                                     ConeH cone);
 
+std::optional<ParamPoint> findVertex(Matrix<MPInt> equations);
+
+std::vector<std::pair<PresburgerRelation, GeneratingFunction>> polytopeGeneratingFunction(PolyhedronH poly);
+
 /// Find a vector that is not orthogonal to any of the given vectors,
 /// i.e., has nonzero dot product with those of the given vectors
 /// that are not null.
diff --git a/mlir/include/mlir/Analysis/Presburger/GeneratingFunction.h b/mlir/include/mlir/Analysis/Presburger/GeneratingFunction.h
index c38eab6efd0fc16..35f15319b4568eb 100644
--- a/mlir/include/mlir/Analysis/Presburger/GeneratingFunction.h
+++ b/mlir/include/mlir/Analysis/Presburger/GeneratingFunction.h
@@ -72,7 +72,7 @@ class GeneratingFunction {
     return denominators;
   }
 
-  GeneratingFunction operator+(GeneratingFunction &gf) const {
+  GeneratingFunction operator+(GeneratingFunction gf) const {
     assert(numParam == gf.getNumParams() &&
            "two generating functions with different numbers of parameters "
            "cannot be added!");
diff --git a/mlir/include/mlir/Analysis/Presburger/IntegerRelation.h b/mlir/include/mlir/Analysis/Presburger/IntegerRelation.h
index c476a022a482725..e7b03c44d33e38c 100644
--- a/mlir/include/mlir/Analysis/Presburger/IntegerRelation.h
+++ b/mlir/include/mlir/Analysis/Presburger/IntegerRelation.h
@@ -711,6 +711,10 @@ class IntegerRelation {
   /// return `this \ set`.
   PresburgerRelation subtract(const PresburgerRelation &set) const;
 
+  void removeTrivialEqualities();
+
+  bool isFullDim();
+
   void print(raw_ostream &os) const;
   void dump() const;
 
diff --git a/mlir/lib/Analysis/Presburger/Barvinok.cpp b/mlir/lib/Analysis/Presburger/Barvinok.cpp
index e0fd0dd8caa4d31..c8cec45a21dc4bb 100644
--- a/mlir/lib/Analysis/Presburger/Barvinok.cpp
+++ b/mlir/lib/Analysis/Presburger/Barvinok.cpp
@@ -10,6 +10,7 @@
 #include "mlir/Analysis/Presburger/Utils.h"
 #include "llvm/ADT/Sequence.h"
 #include <algorithm>
+#include <bitset>
 
 using namespace mlir;
 using namespace presburger;
@@ -147,6 +148,222 @@ GeneratingFunction mlir::presburger::detail::unimodularConeGeneratingFunction(
                             std::vector({denominator}));
 }
 
+std::optional<ParamPoint> mlir::presburger::detail::findVertex(Matrix<MPInt> equations)
+{
+    // `equalities` is a d x (d + p + 1) matrix.
+
+    unsigned r = equations.getNumRows();
+    unsigned c = equations.getNumColumns();
+
+    IntMatrix coeffs(r, r);
+    for (unsigned i = 0; i < r; i++)
+        for (unsigned j = 0; j < r; j++)
+            coeffs(i, j) = equations(i, j);
+    
+    if (coeffs.determinant() == MPInt(0))
+        return std::nullopt;
+
+    Matrix<Fraction> equationsF(r, c);
+    for (unsigned i = 0; i < r; i++)
+        for (unsigned j = 0; j < c; j++)
+            equationsF(i, j) = Fraction(equations(i, j), 1);
+    
+    Fraction a, b;
+    for (unsigned i = 0; i < r; i++)
+    {
+        if (equationsF(i, i) == Fraction(0, 1))
+            for (unsigned j = i+1; j < r; j++)
+                if (equationsF(j, i) != 0)
+                {
+                    equationsF.addToRow(i, equationsF.getRow(j), Fraction(1, 1));
+                    break;
+                }
+        b = equationsF(i, i);
+
+        for (unsigned j = 0; j < r; j++)
+        {
+            if (equationsF(j, i) == 0 || j == i) continue;
+            a = equationsF(j, i);
+            equationsF.addToRow(j, equationsF.getRow(i), - a / b);
+        }
+    }
+
+    for (unsigned i = 0; i < r; i++)
+    {
+        a = equationsF(i, i);
+        for (unsigned j = 0; j < c; j++)
+            equationsF(i, j) = equationsF(i, j) / a;
+    }
+
+    ParamPoint vertex(r, c-r); // d x p+1
+    for (unsigned i = 0; i < r; i++)
+        for (unsigned j = 0; j < c-r; j++)
+            vertex(i, j) = -equationsF(i, r+j);
+    
+    return vertex;
+}
+
+std::vector<std::pair<PresburgerRelation, GeneratingFunction>> mlir::presburger::detail::polytopeGeneratingFunction(PolyhedronH poly)
+{
+    unsigned d = poly.getNumRangeVars(); 
+    unsigned p = poly.getNumSymbolVars();
+    unsigned n = poly.getNumInequalities();
+
+    SmallVector<std::pair<int, ConeH>, 4> unimodCones;
+    GeneratingFunction chamberGf(p, {}, {}, {});
+    std::vector<std::pair<PresburgerRelation, GeneratingFunction>> gf({});
+    ConeH tgtCone = defineHRep(d);
+
+    Matrix<MPInt> subset(d, d+p+1);
+    std::vector<Matrix<MPInt>> subsets; // Stores the inequality subsets corresponding to each vertex.
+    Matrix<Fraction> remaining(n-d, d+p+1);
+
+    std::optional<ParamPoint> vertex;
+    std::vector<ParamPoint> vertices;
+
+    Matrix<Fraction> a2(n-d, d);
+    Matrix<Fraction> b2c2(n-d, p+1);
+
+    Matrix<Fraction> activeRegion(n-d, p+1);
+    Matrix<MPInt> activeRegionNorm(n-d, p+1);
+    MPInt lcmDenoms;
+    IntegerRelation activeRegionRel(PresburgerSpace::getRelationSpace(0, p, 0, 0));
+    // The active region will be defined as activeRegionCoeffs @ p + activeRegionConstant ≥ 0.
+    // The active region is a polyhedron in parameter space.
+    std::vector<PresburgerRelation> activeRegions;
+    
+
+    for (std::bitset<16> indicator(((1ul << d)-1ul) << (n-d));
+        indicator.to_ulong() <= ((1ul << d)-1ul) << (n-d);              // d 1's followed by n-d 0's
+        indicator = std::bitset<16>(indicator.to_ulong() - 1))
+    {
+        if (indicator.count() != d)
+            continue;
+
+        subset = Matrix<MPInt>(d, d+p+1);
+        remaining = Matrix<Fraction>(n-d, d+p+1);
+        unsigned j1 = 0, j2 = 0;
+        for (unsigned i = 0; i < n; i++)
+            if (indicator.test(i))
+                subset.setRow(j1++, poly.getInequality(i));
+                // [A1 | B1 | c1]
+            else
+            {
+                for (unsigned k = 0; k < d; k++)
+                    a2(j2, k) = Fraction(poly.atIneq(i, k), 1);
+                for (unsigned k = d; k < d+p+1; k++)
+                    b2c2(j2, k-d) = Fraction(poly.atIneq(i, k), 1);
+                j2++;
+                // [A2 | B2 | c2]
+            }
+
+        vertex = findVertex(subset); // d x (p+1)
+
+        if (vertex == std::nullopt) continue;
+        vertices.push_back(*vertex);
+        subsets.push_back(subset);
+
+        // Region is given by (A2 @ X + B2) p + (A2 @ y + c2) ≥ 0
+        // This is equivt to A2 @ [X | y] + [B2 | c2]
+        // We premultiply [X | y] with each row of A2 and add each row of [B2 | c2].
+        for (unsigned i = 0; i < n-d; i++)
+        {
+            activeRegion.setRow(i, (*vertex).preMultiplyWithRow(a2.getRow(i)));
+            activeRegion.addToRow(i, b2c2.getRow(i), Fraction(1, 1));
+        }
+
+        activeRegionNorm = Matrix<MPInt>(n-d, p+1);
+        activeRegionRel = IntegerRelation(PresburgerSpace::getRelationSpace(0, p, 0, 0));
+        lcmDenoms = 1;
+        for (unsigned i = 0; i < n-d; i++)
+        {
+            for (unsigned j = 0; j < p+1; j++)
+                lcmDenoms = lcm(lcmDenoms, activeRegion(i, j).den);
+            for (unsigned j = 0; j < p+1; j++)
+                activeRegionNorm(i, j) = (activeRegion(i, j) * Fraction(lcmDenoms, 1)).getAsInteger();
+
+            activeRegionRel.addInequality(activeRegionNorm.getRow(i));
+        }
+
+        activeRegions.push_back(PresburgerRelation(activeRegionRel));
+    }
+
+    // Clauss-Loechner chamber decomposition
+    std::vector<std::pair<PresburgerRelation, std::vector<unsigned>>> chambers = 
+        {std::make_pair(activeRegions[0], std::vector({0u}))};
+    std::vector<std::pair<PresburgerRelation, std::vector<unsigned>>> newChambers;
+    for (unsigned j = 1; j < vertices.size(); j++)
+    {
+        newChambers.clear();
+        PresburgerRelation r_j = activeRegions[j];
+        ParamPoint v_j = vertices[j];
+        for (unsigned i = 0; i < chambers.size(); i++)
+        {
+            auto [r_i, v_i] = chambers[i];
+
+            PresburgerRelation intersection = r_i.intersect(r_j);
+            bool isFullDim = false;
+            for (auto disjunct : intersection.getAllDisjuncts())
+                if (disjunct.isFullDim())
+                {
+                    isFullDim = true;
+                    break;
+                }
+            isFullDim = (p == 0) || isFullDim;
+            if (!isFullDim) newChambers.push_back(chambers[i]);
+            else
+            {
+                PresburgerRelation subtraction = r_i.subtract(r_j);
+                newChambers.push_back(std::make_pair(subtraction, v_i));
+
+                v_i.push_back(j);
+                newChambers.push_back(std::make_pair(intersection, v_i));
+            }
+
+        }
+        for (auto chamber : newChambers)
+            r_j = r_j.subtract(chamber.first);
+
+        newChambers.push_back(std::make_pair(r_j, std::vector({j})));
+
+        chambers.clear();
+        for (auto chamber : newChambers)
+        {
+            bool empty = true;
+            for (auto disjunct : chamber.first.getAllDisjuncts())
+                if (!disjunct.isEmpty())
+                {
+                    empty = false;
+                    break;
+                }
+            if (!empty)
+                chambers.push_back(chamber);
+        }
+    }
+
+    SmallVector<MPInt> ineq(d+1);
+    for (auto chamber : chambers)
+    {
+        chamberGf = GeneratingFunction(p, {}, {}, {});
+        for (unsigned i : chamber.second)
+        {
+            tgtCone = defineHRep(d);
+            for (unsigned j = 0; j < d; j++)
+            {
+                for (unsigned k = 0; k < d; k++)
+                    ineq[k] = subsets[i](j, k);
+                ineq[d] = subsets[i](j, d+p);
+                tgtCone.addInequality(ineq);
+            }
+            unimodCones = {std::make_pair(1, tgtCone)};
+            for (auto signedCone : unimodCones)
+                chamberGf = chamberGf + unimodularConeGeneratingFunction(vertices[i], signedCone.first, signedCone.second);
+        }
+        gf.push_back(std::make_pair(chamber.first, chamberGf));
+    }
+    return gf;
+}
+
 /// We use an iterative procedure to find a vector not orthogonal
 /// to a given set, ignoring the null vectors.
 /// Let the inputs be {x_1, ..., x_k}, all vectors of length n.
diff --git a/mlir/lib/Analysis/Presburger/IntegerRelation.cpp b/mlir/lib/Analysis/Presburger/IntegerRelation.cpp
index 7d2a63d17676f57..578b95571fff242 100644
--- a/mlir/lib/Analysis/Presburger/IntegerRelation.cpp
+++ b/mlir/lib/Analysis/Presburger/IntegerRelation.cpp
@@ -2498,6 +2498,42 @@ void IntegerRelation::printSpace(raw_ostream &os) const {
   os << getNumConstraints() << " constraints\n";
 }
 
+void IntegerRelation::removeTrivialEqualities() {
+    bool flag;
+    for (unsigned i = 0; i < getNumEqualities(); i++)
+    {
+        flag = true;
+        for (unsigned j = 0; j < getNumVars()+1; j++)
+            if (atEq(i, j) != 0)
+                flag = false;
+        if (flag)
+            removeEquality(i);
+    }
+}
+
+bool IntegerRelation::isFullDim() {
+    removeTrivialEqualities(); // to implement: remove equalities that are `0 = 0`
+
+    if (getNumEqualities() > 0)
+      return true;
+
+    Simplex simplex(*this);
+    for (unsigned i = 0; i < getNumInequalities(); i++)
+    {
+      auto ineq = inequalities.getRow(i);
+      auto upOpt = simplex.computeOptimum(Simplex::Direction::Up, ineq);
+      auto downOpt = simplex.computeOptimum(Simplex::Direction::Down, ineq);
+      if (upOpt.getKind() == OptimumKind::Unbounded ||
+          downOpt.getKind() == OptimumKind::Unbounded)
+            continue;
+      if (upOpt.getKind() == OptimumKind::Bounded &&
+          downOpt.getKind() == OptimumKind::Bounded &&
+          *upOpt == *downOpt) // ineq == 0 holds for this
+        return false;
+    }
+    return true;
+  }
+
 void IntegerRelation::print(raw_ostream &os) const {
   assert(hasConsistentState());
   printSpace(os);

>From 3a147d0ee7542720acea7349d15968a1414cf017 Mon Sep 17 00:00:00 2001
From: Abhinav271828 <abhinav.m at research.iiit.ac.in>
Date: Mon, 22 Jan 2024 14:57:29 +0530
Subject: [PATCH 02/26] Formatting

---
 .../mlir/Analysis/Presburger/Barvinok.h       |   3 +-
 mlir/lib/Analysis/Presburger/Barvinok.cpp     | 390 +++++++++---------
 .../Analysis/Presburger/IntegerRelation.cpp   |  54 ++-
 3 files changed, 219 insertions(+), 228 deletions(-)

diff --git a/mlir/include/mlir/Analysis/Presburger/Barvinok.h b/mlir/include/mlir/Analysis/Presburger/Barvinok.h
index 6feb7e8f09233f6..1fafa1e709d5e32 100644
--- a/mlir/include/mlir/Analysis/Presburger/Barvinok.h
+++ b/mlir/include/mlir/Analysis/Presburger/Barvinok.h
@@ -87,7 +87,8 @@ GeneratingFunction unimodularConeGeneratingFunction(ParamPoint vertex, int sign,
 
 std::optional<ParamPoint> findVertex(Matrix<MPInt> equations);
 
-std::vector<std::pair<PresburgerRelation, GeneratingFunction>> polytopeGeneratingFunction(PolyhedronH poly);
+std::vector<std::pair<PresburgerRelation, GeneratingFunction>>
+polytopeGeneratingFunction(PolyhedronH poly);
 
 /// Find a vector that is not orthogonal to any of the given vectors,
 /// i.e., has nonzero dot product with those of the given vectors
diff --git a/mlir/lib/Analysis/Presburger/Barvinok.cpp b/mlir/lib/Analysis/Presburger/Barvinok.cpp
index c8cec45a21dc4bb..f9f75a8a0dcef3b 100644
--- a/mlir/lib/Analysis/Presburger/Barvinok.cpp
+++ b/mlir/lib/Analysis/Presburger/Barvinok.cpp
@@ -148,220 +148,212 @@ GeneratingFunction mlir::presburger::detail::unimodularConeGeneratingFunction(
                             std::vector({denominator}));
 }
 
-std::optional<ParamPoint> mlir::presburger::detail::findVertex(Matrix<MPInt> equations)
-{
-    // `equalities` is a d x (d + p + 1) matrix.
-
-    unsigned r = equations.getNumRows();
-    unsigned c = equations.getNumColumns();
-
-    IntMatrix coeffs(r, r);
-    for (unsigned i = 0; i < r; i++)
-        for (unsigned j = 0; j < r; j++)
-            coeffs(i, j) = equations(i, j);
-    
-    if (coeffs.determinant() == MPInt(0))
-        return std::nullopt;
-
-    Matrix<Fraction> equationsF(r, c);
-    for (unsigned i = 0; i < r; i++)
-        for (unsigned j = 0; j < c; j++)
-            equationsF(i, j) = Fraction(equations(i, j), 1);
-    
-    Fraction a, b;
-    for (unsigned i = 0; i < r; i++)
-    {
-        if (equationsF(i, i) == Fraction(0, 1))
-            for (unsigned j = i+1; j < r; j++)
-                if (equationsF(j, i) != 0)
-                {
-                    equationsF.addToRow(i, equationsF.getRow(j), Fraction(1, 1));
-                    break;
-                }
-        b = equationsF(i, i);
-
-        for (unsigned j = 0; j < r; j++)
-        {
-            if (equationsF(j, i) == 0 || j == i) continue;
-            a = equationsF(j, i);
-            equationsF.addToRow(j, equationsF.getRow(i), - a / b);
+std::optional<ParamPoint>
+mlir::presburger::detail::findVertex(Matrix<MPInt> equations) {
+  // `equalities` is a d x (d + p + 1) matrix.
+
+  unsigned r = equations.getNumRows();
+  unsigned c = equations.getNumColumns();
+
+  IntMatrix coeffs(r, r);
+  for (unsigned i = 0; i < r; i++)
+    for (unsigned j = 0; j < r; j++)
+      coeffs(i, j) = equations(i, j);
+
+  if (coeffs.determinant() == MPInt(0))
+    return std::nullopt;
+
+  Matrix<Fraction> equationsF(r, c);
+  for (unsigned i = 0; i < r; i++)
+    for (unsigned j = 0; j < c; j++)
+      equationsF(i, j) = Fraction(equations(i, j), 1);
+
+  Fraction a, b;
+  for (unsigned i = 0; i < r; i++) {
+    if (equationsF(i, i) == Fraction(0, 1))
+      for (unsigned j = i + 1; j < r; j++)
+        if (equationsF(j, i) != 0) {
+          equationsF.addToRow(i, equationsF.getRow(j), Fraction(1, 1));
+          break;
         }
-    }
+    b = equationsF(i, i);
 
-    for (unsigned i = 0; i < r; i++)
-    {
-        a = equationsF(i, i);
-        for (unsigned j = 0; j < c; j++)
-            equationsF(i, j) = equationsF(i, j) / a;
+    for (unsigned j = 0; j < r; j++) {
+      if (equationsF(j, i) == 0 || j == i)
+        continue;
+      a = equationsF(j, i);
+      equationsF.addToRow(j, equationsF.getRow(i), -a / b);
     }
+  }
 
-    ParamPoint vertex(r, c-r); // d x p+1
-    for (unsigned i = 0; i < r; i++)
-        for (unsigned j = 0; j < c-r; j++)
-            vertex(i, j) = -equationsF(i, r+j);
-    
-    return vertex;
-}
+  for (unsigned i = 0; i < r; i++) {
+    a = equationsF(i, i);
+    for (unsigned j = 0; j < c; j++)
+      equationsF(i, j) = equationsF(i, j) / a;
+  }
 
-std::vector<std::pair<PresburgerRelation, GeneratingFunction>> mlir::presburger::detail::polytopeGeneratingFunction(PolyhedronH poly)
-{
-    unsigned d = poly.getNumRangeVars(); 
-    unsigned p = poly.getNumSymbolVars();
-    unsigned n = poly.getNumInequalities();
-
-    SmallVector<std::pair<int, ConeH>, 4> unimodCones;
-    GeneratingFunction chamberGf(p, {}, {}, {});
-    std::vector<std::pair<PresburgerRelation, GeneratingFunction>> gf({});
-    ConeH tgtCone = defineHRep(d);
-
-    Matrix<MPInt> subset(d, d+p+1);
-    std::vector<Matrix<MPInt>> subsets; // Stores the inequality subsets corresponding to each vertex.
-    Matrix<Fraction> remaining(n-d, d+p+1);
-
-    std::optional<ParamPoint> vertex;
-    std::vector<ParamPoint> vertices;
-
-    Matrix<Fraction> a2(n-d, d);
-    Matrix<Fraction> b2c2(n-d, p+1);
-
-    Matrix<Fraction> activeRegion(n-d, p+1);
-    Matrix<MPInt> activeRegionNorm(n-d, p+1);
-    MPInt lcmDenoms;
-    IntegerRelation activeRegionRel(PresburgerSpace::getRelationSpace(0, p, 0, 0));
-    // The active region will be defined as activeRegionCoeffs @ p + activeRegionConstant ≥ 0.
-    // The active region is a polyhedron in parameter space.
-    std::vector<PresburgerRelation> activeRegions;
-    
-
-    for (std::bitset<16> indicator(((1ul << d)-1ul) << (n-d));
-        indicator.to_ulong() <= ((1ul << d)-1ul) << (n-d);              // d 1's followed by n-d 0's
-        indicator = std::bitset<16>(indicator.to_ulong() - 1))
-    {
-        if (indicator.count() != d)
-            continue;
-
-        subset = Matrix<MPInt>(d, d+p+1);
-        remaining = Matrix<Fraction>(n-d, d+p+1);
-        unsigned j1 = 0, j2 = 0;
-        for (unsigned i = 0; i < n; i++)
-            if (indicator.test(i))
-                subset.setRow(j1++, poly.getInequality(i));
-                // [A1 | B1 | c1]
-            else
-            {
-                for (unsigned k = 0; k < d; k++)
-                    a2(j2, k) = Fraction(poly.atIneq(i, k), 1);
-                for (unsigned k = d; k < d+p+1; k++)
-                    b2c2(j2, k-d) = Fraction(poly.atIneq(i, k), 1);
-                j2++;
-                // [A2 | B2 | c2]
-            }
-
-        vertex = findVertex(subset); // d x (p+1)
-
-        if (vertex == std::nullopt) continue;
-        vertices.push_back(*vertex);
-        subsets.push_back(subset);
-
-        // Region is given by (A2 @ X + B2) p + (A2 @ y + c2) ≥ 0
-        // This is equivt to A2 @ [X | y] + [B2 | c2]
-        // We premultiply [X | y] with each row of A2 and add each row of [B2 | c2].
-        for (unsigned i = 0; i < n-d; i++)
-        {
-            activeRegion.setRow(i, (*vertex).preMultiplyWithRow(a2.getRow(i)));
-            activeRegion.addToRow(i, b2c2.getRow(i), Fraction(1, 1));
-        }
+  ParamPoint vertex(r, c - r); // d x p+1
+  for (unsigned i = 0; i < r; i++)
+    for (unsigned j = 0; j < c - r; j++)
+      vertex(i, j) = -equationsF(i, r + j);
 
-        activeRegionNorm = Matrix<MPInt>(n-d, p+1);
-        activeRegionRel = IntegerRelation(PresburgerSpace::getRelationSpace(0, p, 0, 0));
-        lcmDenoms = 1;
-        for (unsigned i = 0; i < n-d; i++)
-        {
-            for (unsigned j = 0; j < p+1; j++)
-                lcmDenoms = lcm(lcmDenoms, activeRegion(i, j).den);
-            for (unsigned j = 0; j < p+1; j++)
-                activeRegionNorm(i, j) = (activeRegion(i, j) * Fraction(lcmDenoms, 1)).getAsInteger();
-
-            activeRegionRel.addInequality(activeRegionNorm.getRow(i));
-        }
+  return vertex;
+}
 
-        activeRegions.push_back(PresburgerRelation(activeRegionRel));
+std::vector<std::pair<PresburgerRelation, GeneratingFunction>>
+mlir::presburger::detail::polytopeGeneratingFunction(PolyhedronH poly) {
+  unsigned d = poly.getNumRangeVars();
+  unsigned p = poly.getNumSymbolVars();
+  unsigned n = poly.getNumInequalities();
+
+  SmallVector<std::pair<int, ConeH>, 4> unimodCones;
+  GeneratingFunction chamberGf(p, {}, {}, {});
+  std::vector<std::pair<PresburgerRelation, GeneratingFunction>> gf({});
+  ConeH tgtCone = defineHRep(d);
+
+  Matrix<MPInt> subset(d, d + p + 1);
+  std::vector<Matrix<MPInt>>
+      subsets; // Stores the inequality subsets corresponding to each vertex.
+  Matrix<Fraction> remaining(n - d, d + p + 1);
+
+  std::optional<ParamPoint> vertex;
+  std::vector<ParamPoint> vertices;
+
+  Matrix<Fraction> a2(n - d, d);
+  Matrix<Fraction> b2c2(n - d, p + 1);
+
+  Matrix<Fraction> activeRegion(n - d, p + 1);
+  Matrix<MPInt> activeRegionNorm(n - d, p + 1);
+  MPInt lcmDenoms;
+  IntegerRelation activeRegionRel(
+      PresburgerSpace::getRelationSpace(0, p, 0, 0));
+  // The active region will be defined as activeRegionCoeffs @ p +
+  // activeRegionConstant ≥ 0. The active region is a polyhedron in parameter
+  // space.
+  std::vector<PresburgerRelation> activeRegions;
+
+  for (std::bitset<16> indicator(((1ul << d) - 1ul) << (n - d));
+       indicator.to_ulong() <= ((1ul << d) - 1ul)
+                                   << (n - d); // d 1's followed by n-d 0's
+       indicator = std::bitset<16>(indicator.to_ulong() - 1)) {
+    if (indicator.count() != d)
+      continue;
+
+    subset = Matrix<MPInt>(d, d + p + 1);
+    remaining = Matrix<Fraction>(n - d, d + p + 1);
+    unsigned j1 = 0, j2 = 0;
+    for (unsigned i = 0; i < n; i++)
+      if (indicator.test(i))
+        subset.setRow(j1++, poly.getInequality(i));
+      // [A1 | B1 | c1]
+      else {
+        for (unsigned k = 0; k < d; k++)
+          a2(j2, k) = Fraction(poly.atIneq(i, k), 1);
+        for (unsigned k = d; k < d + p + 1; k++)
+          b2c2(j2, k - d) = Fraction(poly.atIneq(i, k), 1);
+        j2++;
+        // [A2 | B2 | c2]
+      }
+
+    vertex = findVertex(subset); // d x (p+1)
+
+    if (vertex == std::nullopt)
+      continue;
+    vertices.push_back(*vertex);
+    subsets.push_back(subset);
+
+    // Region is given by (A2 @ X + B2) p + (A2 @ y + c2) ≥ 0
+    // This is equivt to A2 @ [X | y] + [B2 | c2]
+    // We premultiply [X | y] with each row of A2 and add each row of [B2 | c2].
+    for (unsigned i = 0; i < n - d; i++) {
+      activeRegion.setRow(i, (*vertex).preMultiplyWithRow(a2.getRow(i)));
+      activeRegion.addToRow(i, b2c2.getRow(i), Fraction(1, 1));
     }
 
-    // Clauss-Loechner chamber decomposition
-    std::vector<std::pair<PresburgerRelation, std::vector<unsigned>>> chambers = 
-        {std::make_pair(activeRegions[0], std::vector({0u}))};
-    std::vector<std::pair<PresburgerRelation, std::vector<unsigned>>> newChambers;
-    for (unsigned j = 1; j < vertices.size(); j++)
-    {
-        newChambers.clear();
-        PresburgerRelation r_j = activeRegions[j];
-        ParamPoint v_j = vertices[j];
-        for (unsigned i = 0; i < chambers.size(); i++)
-        {
-            auto [r_i, v_i] = chambers[i];
-
-            PresburgerRelation intersection = r_i.intersect(r_j);
-            bool isFullDim = false;
-            for (auto disjunct : intersection.getAllDisjuncts())
-                if (disjunct.isFullDim())
-                {
-                    isFullDim = true;
-                    break;
-                }
-            isFullDim = (p == 0) || isFullDim;
-            if (!isFullDim) newChambers.push_back(chambers[i]);
-            else
-            {
-                PresburgerRelation subtraction = r_i.subtract(r_j);
-                newChambers.push_back(std::make_pair(subtraction, v_i));
-
-                v_i.push_back(j);
-                newChambers.push_back(std::make_pair(intersection, v_i));
-            }
+    activeRegionNorm = Matrix<MPInt>(n - d, p + 1);
+    activeRegionRel =
+        IntegerRelation(PresburgerSpace::getRelationSpace(0, p, 0, 0));
+    lcmDenoms = 1;
+    for (unsigned i = 0; i < n - d; i++) {
+      for (unsigned j = 0; j < p + 1; j++)
+        lcmDenoms = lcm(lcmDenoms, activeRegion(i, j).den);
+      for (unsigned j = 0; j < p + 1; j++)
+        activeRegionNorm(i, j) =
+            (activeRegion(i, j) * Fraction(lcmDenoms, 1)).getAsInteger();
+
+      activeRegionRel.addInequality(activeRegionNorm.getRow(i));
+    }
+
+    activeRegions.push_back(PresburgerRelation(activeRegionRel));
+  }
 
+  // Clauss-Loechner chamber decomposition
+  std::vector<std::pair<PresburgerRelation, std::vector<unsigned>>> chambers = {
+      std::make_pair(activeRegions[0], std::vector({0u}))};
+  std::vector<std::pair<PresburgerRelation, std::vector<unsigned>>> newChambers;
+  for (unsigned j = 1; j < vertices.size(); j++) {
+    newChambers.clear();
+    PresburgerRelation r_j = activeRegions[j];
+    ParamPoint v_j = vertices[j];
+    for (unsigned i = 0; i < chambers.size(); i++) {
+      auto [r_i, v_i] = chambers[i];
+
+      PresburgerRelation intersection = r_i.intersect(r_j);
+      bool isFullDim = false;
+      for (auto disjunct : intersection.getAllDisjuncts())
+        if (disjunct.isFullDim()) {
+          isFullDim = true;
+          break;
         }
-        for (auto chamber : newChambers)
-            r_j = r_j.subtract(chamber.first);
-
-        newChambers.push_back(std::make_pair(r_j, std::vector({j})));
-
-        chambers.clear();
-        for (auto chamber : newChambers)
-        {
-            bool empty = true;
-            for (auto disjunct : chamber.first.getAllDisjuncts())
-                if (!disjunct.isEmpty())
-                {
-                    empty = false;
-                    break;
-                }
-            if (!empty)
-                chambers.push_back(chamber);
+      isFullDim = (p == 0) || isFullDim;
+      if (!isFullDim)
+        newChambers.push_back(chambers[i]);
+      else {
+        PresburgerRelation subtraction = r_i.subtract(r_j);
+        newChambers.push_back(std::make_pair(subtraction, v_i));
+
+        v_i.push_back(j);
+        newChambers.push_back(std::make_pair(intersection, v_i));
+      }
+    }
+    for (auto chamber : newChambers)
+      r_j = r_j.subtract(chamber.first);
+
+    newChambers.push_back(std::make_pair(r_j, std::vector({j})));
+
+    chambers.clear();
+    for (auto chamber : newChambers) {
+      bool empty = true;
+      for (auto disjunct : chamber.first.getAllDisjuncts())
+        if (!disjunct.isEmpty()) {
+          empty = false;
+          break;
         }
+      if (!empty)
+        chambers.push_back(chamber);
     }
+  }
 
-    SmallVector<MPInt> ineq(d+1);
-    for (auto chamber : chambers)
-    {
-        chamberGf = GeneratingFunction(p, {}, {}, {});
-        for (unsigned i : chamber.second)
-        {
-            tgtCone = defineHRep(d);
-            for (unsigned j = 0; j < d; j++)
-            {
-                for (unsigned k = 0; k < d; k++)
-                    ineq[k] = subsets[i](j, k);
-                ineq[d] = subsets[i](j, d+p);
-                tgtCone.addInequality(ineq);
-            }
-            unimodCones = {std::make_pair(1, tgtCone)};
-            for (auto signedCone : unimodCones)
-                chamberGf = chamberGf + unimodularConeGeneratingFunction(vertices[i], signedCone.first, signedCone.second);
-        }
-        gf.push_back(std::make_pair(chamber.first, chamberGf));
+  SmallVector<MPInt> ineq(d + 1);
+  for (auto chamber : chambers) {
+    chamberGf = GeneratingFunction(p, {}, {}, {});
+    for (unsigned i : chamber.second) {
+      tgtCone = defineHRep(d);
+      for (unsigned j = 0; j < d; j++) {
+        for (unsigned k = 0; k < d; k++)
+          ineq[k] = subsets[i](j, k);
+        ineq[d] = subsets[i](j, d + p);
+        tgtCone.addInequality(ineq);
+      }
+      unimodCones = {std::make_pair(1, tgtCone)};
+      for (auto signedCone : unimodCones)
+        chamberGf =
+            chamberGf + unimodularConeGeneratingFunction(
+                            vertices[i], signedCone.first, signedCone.second);
     }
-    return gf;
+    gf.push_back(std::make_pair(chamber.first, chamberGf));
+  }
+  return gf;
 }
 
 /// We use an iterative procedure to find a vector not orthogonal
diff --git a/mlir/lib/Analysis/Presburger/IntegerRelation.cpp b/mlir/lib/Analysis/Presburger/IntegerRelation.cpp
index 578b95571fff242..edb95b4072e1d85 100644
--- a/mlir/lib/Analysis/Presburger/IntegerRelation.cpp
+++ b/mlir/lib/Analysis/Presburger/IntegerRelation.cpp
@@ -2499,40 +2499,38 @@ void IntegerRelation::printSpace(raw_ostream &os) const {
 }
 
 void IntegerRelation::removeTrivialEqualities() {
-    bool flag;
-    for (unsigned i = 0; i < getNumEqualities(); i++)
-    {
-        flag = true;
-        for (unsigned j = 0; j < getNumVars()+1; j++)
-            if (atEq(i, j) != 0)
-                flag = false;
-        if (flag)
-            removeEquality(i);
-    }
+  bool flag;
+  for (unsigned i = 0; i < getNumEqualities(); i++) {
+    flag = true;
+    for (unsigned j = 0; j < getNumVars() + 1; j++)
+      if (atEq(i, j) != 0)
+        flag = false;
+    if (flag)
+      removeEquality(i);
+  }
 }
 
 bool IntegerRelation::isFullDim() {
-    removeTrivialEqualities(); // to implement: remove equalities that are `0 = 0`
-
-    if (getNumEqualities() > 0)
-      return true;
+  removeTrivialEqualities(); // to implement: remove equalities that are `0 = 0`
 
-    Simplex simplex(*this);
-    for (unsigned i = 0; i < getNumInequalities(); i++)
-    {
-      auto ineq = inequalities.getRow(i);
-      auto upOpt = simplex.computeOptimum(Simplex::Direction::Up, ineq);
-      auto downOpt = simplex.computeOptimum(Simplex::Direction::Down, ineq);
-      if (upOpt.getKind() == OptimumKind::Unbounded ||
-          downOpt.getKind() == OptimumKind::Unbounded)
-            continue;
-      if (upOpt.getKind() == OptimumKind::Bounded &&
-          downOpt.getKind() == OptimumKind::Bounded &&
-          *upOpt == *downOpt) // ineq == 0 holds for this
-        return false;
-    }
+  if (getNumEqualities() > 0)
     return true;
+
+  Simplex simplex(*this);
+  for (unsigned i = 0; i < getNumInequalities(); i++) {
+    auto ineq = inequalities.getRow(i);
+    auto upOpt = simplex.computeOptimum(Simplex::Direction::Up, ineq);
+    auto downOpt = simplex.computeOptimum(Simplex::Direction::Down, ineq);
+    if (upOpt.getKind() == OptimumKind::Unbounded ||
+        downOpt.getKind() == OptimumKind::Unbounded)
+      continue;
+    if (upOpt.getKind() == OptimumKind::Bounded &&
+        downOpt.getKind() == OptimumKind::Bounded &&
+        *upOpt == *downOpt) // ineq == 0 holds for this
+      return false;
   }
+  return true;
+}
 
 void IntegerRelation::print(raw_ostream &os) const {
   assert(hasConsistentState());

>From 294dc0ae9147463248b9b40f04d4a2529ff09ed4 Mon Sep 17 00:00:00 2001
From: Abhinav271828 <abhinav.m at research.iiit.ac.in>
Date: Mon, 22 Jan 2024 15:07:38 +0530
Subject: [PATCH 03/26] Use intmatrix and fracmatrix

---
 .../mlir/Analysis/Presburger/Barvinok.h       |  2 +-
 mlir/lib/Analysis/Presburger/Barvinok.cpp     | 24 +++++++++----------
 2 files changed, 13 insertions(+), 13 deletions(-)

diff --git a/mlir/include/mlir/Analysis/Presburger/Barvinok.h b/mlir/include/mlir/Analysis/Presburger/Barvinok.h
index 1fafa1e709d5e32..320848348552c93 100644
--- a/mlir/include/mlir/Analysis/Presburger/Barvinok.h
+++ b/mlir/include/mlir/Analysis/Presburger/Barvinok.h
@@ -85,7 +85,7 @@ ConeH getDual(ConeV cone);
 GeneratingFunction unimodularConeGeneratingFunction(ParamPoint vertex, int sign,
                                                     ConeH cone);
 
-std::optional<ParamPoint> findVertex(Matrix<MPInt> equations);
+std::optional<ParamPoint> findVertex(IntMatrix equations);
 
 std::vector<std::pair<PresburgerRelation, GeneratingFunction>>
 polytopeGeneratingFunction(PolyhedronH poly);
diff --git a/mlir/lib/Analysis/Presburger/Barvinok.cpp b/mlir/lib/Analysis/Presburger/Barvinok.cpp
index f9f75a8a0dcef3b..3753f715c7d925b 100644
--- a/mlir/lib/Analysis/Presburger/Barvinok.cpp
+++ b/mlir/lib/Analysis/Presburger/Barvinok.cpp
@@ -149,7 +149,7 @@ GeneratingFunction mlir::presburger::detail::unimodularConeGeneratingFunction(
 }
 
 std::optional<ParamPoint>
-mlir::presburger::detail::findVertex(Matrix<MPInt> equations) {
+mlir::presburger::detail::findVertex(IntMatrix equations) {
   // `equalities` is a d x (d + p + 1) matrix.
 
   unsigned r = equations.getNumRows();
@@ -163,7 +163,7 @@ mlir::presburger::detail::findVertex(Matrix<MPInt> equations) {
   if (coeffs.determinant() == MPInt(0))
     return std::nullopt;
 
-  Matrix<Fraction> equationsF(r, c);
+  FracMatrix equationsF(r, c);
   for (unsigned i = 0; i < r; i++)
     for (unsigned j = 0; j < c; j++)
       equationsF(i, j) = Fraction(equations(i, j), 1);
@@ -211,19 +211,19 @@ mlir::presburger::detail::polytopeGeneratingFunction(PolyhedronH poly) {
   std::vector<std::pair<PresburgerRelation, GeneratingFunction>> gf({});
   ConeH tgtCone = defineHRep(d);
 
-  Matrix<MPInt> subset(d, d + p + 1);
-  std::vector<Matrix<MPInt>>
+  IntMatrix subset(d, d + p + 1);
+  std::vector<IntMatrix>
       subsets; // Stores the inequality subsets corresponding to each vertex.
-  Matrix<Fraction> remaining(n - d, d + p + 1);
+  FracMatrix remaining(n - d, d + p + 1);
 
   std::optional<ParamPoint> vertex;
   std::vector<ParamPoint> vertices;
 
-  Matrix<Fraction> a2(n - d, d);
-  Matrix<Fraction> b2c2(n - d, p + 1);
+  FracMatrix a2(n - d, d);
+  FracMatrix b2c2(n - d, p + 1);
 
-  Matrix<Fraction> activeRegion(n - d, p + 1);
-  Matrix<MPInt> activeRegionNorm(n - d, p + 1);
+  FracMatrix activeRegion(n - d, p + 1);
+  IntMatrix activeRegionNorm(n - d, p + 1);
   MPInt lcmDenoms;
   IntegerRelation activeRegionRel(
       PresburgerSpace::getRelationSpace(0, p, 0, 0));
@@ -239,8 +239,8 @@ mlir::presburger::detail::polytopeGeneratingFunction(PolyhedronH poly) {
     if (indicator.count() != d)
       continue;
 
-    subset = Matrix<MPInt>(d, d + p + 1);
-    remaining = Matrix<Fraction>(n - d, d + p + 1);
+    subset = IntMatrix(d, d + p + 1);
+    remaining = FracMatrix(n - d, d + p + 1);
     unsigned j1 = 0, j2 = 0;
     for (unsigned i = 0; i < n; i++)
       if (indicator.test(i))
@@ -270,7 +270,7 @@ mlir::presburger::detail::polytopeGeneratingFunction(PolyhedronH poly) {
       activeRegion.addToRow(i, b2c2.getRow(i), Fraction(1, 1));
     }
 
-    activeRegionNorm = Matrix<MPInt>(n - d, p + 1);
+    activeRegionNorm = IntMatrix(n - d, p + 1);
     activeRegionRel =
         IntegerRelation(PresburgerSpace::getRelationSpace(0, p, 0, 0));
     lcmDenoms = 1;

>From 35e9b7665fcb8b5445d39aa98274e3e63904e368 Mon Sep 17 00:00:00 2001
From: Abhinav271828 <abhinav.m at research.iiit.ac.in>
Date: Mon, 22 Jan 2024 19:11:47 +0530
Subject: [PATCH 04/26] Doc

---
 .../mlir/Analysis/Presburger/Barvinok.h       |  4 ++
 mlir/lib/Analysis/Presburger/Barvinok.cpp     | 71 ++++++++++++-------
 2 files changed, 50 insertions(+), 25 deletions(-)

diff --git a/mlir/include/mlir/Analysis/Presburger/Barvinok.h b/mlir/include/mlir/Analysis/Presburger/Barvinok.h
index 320848348552c93..9c940231d5e1471 100644
--- a/mlir/include/mlir/Analysis/Presburger/Barvinok.h
+++ b/mlir/include/mlir/Analysis/Presburger/Barvinok.h
@@ -85,6 +85,10 @@ ConeH getDual(ConeV cone);
 GeneratingFunction unimodularConeGeneratingFunction(ParamPoint vertex, int sign,
                                                     ConeH cone);
 
+/// Find the solution of a set of equations that express affine constraints
+/// between a set of variables and a set of parameters. The solution expresses
+/// each variable as an affine function of the parameters.
+/// If there is no solution, return std::null.
 std::optional<ParamPoint> findVertex(IntMatrix equations);
 
 std::vector<std::pair<PresburgerRelation, GeneratingFunction>>
diff --git a/mlir/lib/Analysis/Presburger/Barvinok.cpp b/mlir/lib/Analysis/Presburger/Barvinok.cpp
index 3753f715c7d925b..ff202984bd33624 100644
--- a/mlir/lib/Analysis/Presburger/Barvinok.cpp
+++ b/mlir/lib/Analysis/Presburger/Barvinok.cpp
@@ -148,54 +148,75 @@ GeneratingFunction mlir::presburger::detail::unimodularConeGeneratingFunction(
                             std::vector({denominator}));
 }
 
+/// We use Gaussian elimination to find the solution to a set of d equations
+/// of the form
+/// a_1 x_1 + ... + a_d x_d + b_1 m_1 + ... + b_p m_p + c = 0
+/// where x_i are variables,
+/// m_i are parameters and
+/// a_i, b_i, c are rational coefficients.
+/// The solution expresses each x_i as an affine function of the m_i, and is
+/// therefore represented as a matrix of size d x (p+1).
+/// If there is no solution, we return null.
 std::optional<ParamPoint>
 mlir::presburger::detail::findVertex(IntMatrix equations) {
-  // `equalities` is a d x (d + p + 1) matrix.
+  // equations is a d x (d + p + 1) matrix.
+  // Each row represents an equation.
 
-  unsigned r = equations.getNumRows();
+  unsigned d = equations.getNumRows();
   unsigned c = equations.getNumColumns();
 
-  IntMatrix coeffs(r, r);
-  for (unsigned i = 0; i < r; i++)
-    for (unsigned j = 0; j < r; j++)
+  // First, we check that the system has a solution, and return
+  // null if not.
+  IntMatrix coeffs(d, d);
+  for (unsigned i = 0; i < d; i++)
+    for (unsigned j = 0; j < d; j++)
       coeffs(i, j) = equations(i, j);
 
-  if (coeffs.determinant() == MPInt(0))
+  if (coeffs.determinant() == 0)
     return std::nullopt;
 
-  FracMatrix equationsF(r, c);
-  for (unsigned i = 0; i < r; i++)
-    for (unsigned j = 0; j < c; j++)
-      equationsF(i, j) = Fraction(equations(i, j), 1);
+  // We work with rational numbers.
+  FracMatrix equationsF(equations);
 
-  Fraction a, b;
-  for (unsigned i = 0; i < r; i++) {
-    if (equationsF(i, i) == Fraction(0, 1))
-      for (unsigned j = i + 1; j < r; j++)
+  for (unsigned i = 0; i < d; ++i) {
+    // First ensure that the diagonal element is nonzero, by swapping
+    // it with a nonzero row.
+    if (equationsF(i, i) == 0) {
+      for (unsigned j = i + 1; j < d; ++j) {
         if (equationsF(j, i) != 0) {
-          equationsF.addToRow(i, equationsF.getRow(j), Fraction(1, 1));
+          equationsF.swapRows(j, i);
           break;
         }
-    b = equationsF(i, i);
+      }
+    }
+
+    Fraction b = equationsF(i, i);
 
-    for (unsigned j = 0; j < r; j++) {
+    // Set all elements except the diagonal to zero.
+    for (unsigned j = 0; j < d; ++j) {
       if (equationsF(j, i) == 0 || j == i)
         continue;
-      a = equationsF(j, i);
+      // Set element (j, i) to zero
+      // by subtracting the ith row,
+      // appropriately scaled.
+      Fraction a = equationsF(j, i);
       equationsF.addToRow(j, equationsF.getRow(i), -a / b);
     }
   }
 
-  for (unsigned i = 0; i < r; i++) {
-    a = equationsF(i, i);
-    for (unsigned j = 0; j < c; j++)
+  // Rescale diagonal elements to 1.
+  for (unsigned i = 0; i < d; ++i) {
+    Fraction a = equationsF(i, i);
+    for (unsigned j = 0; j < c; ++j)
       equationsF(i, j) = equationsF(i, j) / a;
   }
 
-  ParamPoint vertex(r, c - r); // d x p+1
-  for (unsigned i = 0; i < r; i++)
-    for (unsigned j = 0; j < c - r; j++)
-      vertex(i, j) = -equationsF(i, r + j);
+  // We copy the last p+1 columns of the matrix as the values of x_i.
+  // We shift the parameter terms to the RHS, and so flip their sign.
+  ParamPoint vertex(d, c - d);
+  for (unsigned i = 0; i < d; ++i)
+    for (unsigned j = 0; j < c - d; ++j)
+      vertex(i, j) = -equationsF(i, d + j);
 
   return vertex;
 }

>From 15e20fe3d2fe845efcb6b3f32165117980ce46dd Mon Sep 17 00:00:00 2001
From: Abhinav271828 <abhinav.m at research.iiit.ac.in>
Date: Mon, 22 Jan 2024 20:19:08 +0530
Subject: [PATCH 05/26] Doc

---
 .../mlir/Analysis/Presburger/Barvinok.h       |   2 +
 mlir/lib/Analysis/Presburger/Barvinok.cpp     | 228 ++++++++++++------
 2 files changed, 157 insertions(+), 73 deletions(-)

diff --git a/mlir/include/mlir/Analysis/Presburger/Barvinok.h b/mlir/include/mlir/Analysis/Presburger/Barvinok.h
index 9c940231d5e1471..a5bc49d79408c6c 100644
--- a/mlir/include/mlir/Analysis/Presburger/Barvinok.h
+++ b/mlir/include/mlir/Analysis/Presburger/Barvinok.h
@@ -91,6 +91,8 @@ GeneratingFunction unimodularConeGeneratingFunction(ParamPoint vertex, int sign,
 /// If there is no solution, return std::null.
 std::optional<ParamPoint> findVertex(IntMatrix equations);
 
+/// Compute the generating function corresponding to a polytope.
+/// All tangent cones of the polytope must be unimodular.
 std::vector<std::pair<PresburgerRelation, GeneratingFunction>>
 polytopeGeneratingFunction(PolyhedronH poly);
 
diff --git a/mlir/lib/Analysis/Presburger/Barvinok.cpp b/mlir/lib/Analysis/Presburger/Barvinok.cpp
index ff202984bd33624..3afc9ad8638ccd2 100644
--- a/mlir/lib/Analysis/Presburger/Barvinok.cpp
+++ b/mlir/lib/Analysis/Presburger/Barvinok.cpp
@@ -162,14 +162,14 @@ mlir::presburger::detail::findVertex(IntMatrix equations) {
   // equations is a d x (d + p + 1) matrix.
   // Each row represents an equation.
 
-  unsigned d = equations.getNumRows();
-  unsigned c = equations.getNumColumns();
+  unsigned numEqs = equations.getNumRows();
+  unsigned numCols = equations.getNumColumns();
 
   // First, we check that the system has a solution, and return
   // null if not.
-  IntMatrix coeffs(d, d);
-  for (unsigned i = 0; i < d; i++)
-    for (unsigned j = 0; j < d; j++)
+  IntMatrix coeffs(numEqs, numEqs);
+  for (unsigned i = 0; i < numEqs; i++)
+    for (unsigned j = 0; j < numEqs; j++)
       coeffs(i, j) = equations(i, j);
 
   if (coeffs.determinant() == 0)
@@ -178,11 +178,11 @@ mlir::presburger::detail::findVertex(IntMatrix equations) {
   // We work with rational numbers.
   FracMatrix equationsF(equations);
 
-  for (unsigned i = 0; i < d; ++i) {
+  for (unsigned i = 0; i < numEqs; ++i) {
     // First ensure that the diagonal element is nonzero, by swapping
     // it with a nonzero row.
     if (equationsF(i, i) == 0) {
-      for (unsigned j = i + 1; j < d; ++j) {
+      for (unsigned j = i + 1; j < numEqs; ++j) {
         if (equationsF(j, i) != 0) {
           equationsF.swapRows(j, i);
           break;
@@ -193,7 +193,7 @@ mlir::presburger::detail::findVertex(IntMatrix equations) {
     Fraction b = equationsF(i, i);
 
     // Set all elements except the diagonal to zero.
-    for (unsigned j = 0; j < d; ++j) {
+    for (unsigned j = 0; j < numEqs; ++j) {
       if (equationsF(j, i) == 0 || j == i)
         continue;
       // Set element (j, i) to zero
@@ -205,102 +205,140 @@ mlir::presburger::detail::findVertex(IntMatrix equations) {
   }
 
   // Rescale diagonal elements to 1.
-  for (unsigned i = 0; i < d; ++i) {
+  for (unsigned i = 0; i < numEqs; ++i) {
     Fraction a = equationsF(i, i);
-    for (unsigned j = 0; j < c; ++j)
+    for (unsigned j = 0; j < numCols; ++j)
       equationsF(i, j) = equationsF(i, j) / a;
   }
 
   // We copy the last p+1 columns of the matrix as the values of x_i.
   // We shift the parameter terms to the RHS, and so flip their sign.
-  ParamPoint vertex(d, c - d);
-  for (unsigned i = 0; i < d; ++i)
-    for (unsigned j = 0; j < c - d; ++j)
-      vertex(i, j) = -equationsF(i, d + j);
+  ParamPoint vertex(numEqs, numCols - numEqs);
+  for (unsigned i = 0; i < numEqs; ++i)
+    for (unsigned j = 0; j < numCols - numEqs; ++j)
+      vertex(i, j) = -equationsF(i, numEqs + j);
 
   return vertex;
 }
 
+/// For a polytope expressed as a set of inequalities, compute the generating
+/// function corresponding to the number of lattice points present. This
+/// algorithm has three main steps:
+/// 1. Enumerate the vertices, by iterating over subsets of inequalities and
+///    checking for solubility.
+/// 2. For each vertex, identify the tangent cone and compute the generating
+///    function corresponding to it. The sum of these GFs is the GF of the
+///    polytope.
+/// 3. [Clauss-Loechner decomposition] Identify the regions in parameter space
+///    (chambers) where each vertex is active, and accordingly compute the
+///    GF of the polytope in each chamber.
+///
+/// Verdoolaege, Sven, et al. "Counting integer points in parametric
+/// polytopes using Barvinok's rational functions." Algorithmica 48 (2007):
+/// 37-66.
 std::vector<std::pair<PresburgerRelation, GeneratingFunction>>
 mlir::presburger::detail::polytopeGeneratingFunction(PolyhedronH poly) {
-  unsigned d = poly.getNumRangeVars();
-  unsigned p = poly.getNumSymbolVars();
-  unsigned n = poly.getNumInequalities();
+  unsigned numVars = poly.getNumRangeVars();
+  unsigned numParams = poly.getNumSymbolVars();
+  unsigned numIneqs = poly.getNumInequalities();
 
-  SmallVector<std::pair<int, ConeH>, 4> unimodCones;
-  GeneratingFunction chamberGf(p, {}, {}, {});
+  // The generating function of the polytope is computed as a set of generating
+  // functions, each one associated with a region in parameter space (chamber).
   std::vector<std::pair<PresburgerRelation, GeneratingFunction>> gf({});
-  ConeH tgtCone = defineHRep(d);
-
-  IntMatrix subset(d, d + p + 1);
-  std::vector<IntMatrix>
-      subsets; // Stores the inequality subsets corresponding to each vertex.
-  FracMatrix remaining(n - d, d + p + 1);
 
-  std::optional<ParamPoint> vertex;
-  std::vector<ParamPoint> vertices;
-
-  FracMatrix a2(n - d, d);
-  FracMatrix b2c2(n - d, p + 1);
-
-  FracMatrix activeRegion(n - d, p + 1);
-  IntMatrix activeRegionNorm(n - d, p + 1);
-  MPInt lcmDenoms;
-  IntegerRelation activeRegionRel(
-      PresburgerSpace::getRelationSpace(0, p, 0, 0));
   // The active region will be defined as activeRegionCoeffs @ p +
   // activeRegionConstant ≥ 0. The active region is a polyhedron in parameter
   // space.
+  FracMatrix activeRegion(numIneqs - numVars, numParams + 1);
+
+  // These vectors store lists of
+  // subsets of inequalities,
+  // the vertices corresponding to them, and
+  // the active regions of the vertices, in order.
+  std::vector<IntMatrix> subsets;
+  std::vector<ParamPoint> vertices;
   std::vector<PresburgerRelation> activeRegions;
 
-  for (std::bitset<16> indicator(((1ul << d) - 1ul) << (n - d));
-       indicator.to_ulong() <= ((1ul << d) - 1ul)
-                                   << (n - d); // d 1's followed by n-d 0's
+  FracMatrix a2(numIneqs - numVars, numVars);
+  FracMatrix b2c2(numIneqs - numVars, numParams + 1);
+
+  // We iterate over all subsets of inequalities with cardinality numVars,
+  // using bitsets up to 2^numIneqs to enumerate.
+  for (std::bitset<16> indicator(((1ul << numVars) - 1ul)
+                                 << (numIneqs - numVars));
+       indicator.to_ulong() <=
+       ((1ul << numVars) - 1ul)
+           << (numIneqs - numVars); // d 1's followed by n-numVars 0's
        indicator = std::bitset<16>(indicator.to_ulong() - 1)) {
-    if (indicator.count() != d)
+
+    if (indicator.count() != numVars)
       continue;
 
-    subset = IntMatrix(d, d + p + 1);
-    remaining = FracMatrix(n - d, d + p + 1);
+    // Collect the inequalities corresponding to the bits which are set.
+    IntMatrix subset(numVars, numVars + numParams + 1);
     unsigned j1 = 0, j2 = 0;
-    for (unsigned i = 0; i < n; i++)
+    for (unsigned i = 0; i < numIneqs; i++)
       if (indicator.test(i))
         subset.setRow(j1++, poly.getInequality(i));
-      // [A1 | B1 | c1]
+
       else {
-        for (unsigned k = 0; k < d; k++)
-          a2(j2, k) = Fraction(poly.atIneq(i, k), 1);
-        for (unsigned k = d; k < d + p + 1; k++)
-          b2c2(j2, k - d) = Fraction(poly.atIneq(i, k), 1);
+        // All other inequalities are stored in a2 and b2c2.
+        // These are column-wise splits of the inequalities;
+        // a2 stores the coefficients of the variables, and
+        // b2c2 stores the coefficients of the parameters and the constant term.
+        for (unsigned k = 0; k < numVars; k++)
+          a2(j2, k) = poly.atIneq(i, k);
+        for (unsigned k = numVars; k < numVars + numParams + 1; k++)
+          b2c2(j2, k - numVars) = poly.atIneq(i, k);
         j2++;
-        // [A2 | B2 | c2]
       }
 
-    vertex = findVertex(subset); // d x (p+1)
+    // Find the vertex, if any, corresponding to the current subset of
+    // inequalities.
+    std::optional<ParamPoint> vertex = findVertex(subset); // d x (p+1)
 
     if (vertex == std::nullopt)
       continue;
+    // If this subset corresponds to a vertex, store it.
     vertices.push_back(*vertex);
     subsets.push_back(subset);
 
-    // Region is given by (A2 @ X + B2) p + (A2 @ y + c2) ≥ 0
-    // This is equivt to A2 @ [X | y] + [B2 | c2]
-    // We premultiply [X | y] with each row of A2 and add each row of [B2 | c2].
-    for (unsigned i = 0; i < n - d; i++) {
+    // Let the current vertex be [X | y], where
+    // X represents the coefficients of the parameters and
+    // y represents the constant term.
+
+    // The region (in parameter space) where this vertex is active is given
+    // by substituting the vertex into the *remaining* inequalities of the
+    // polytope (those which were not collected into `subset`), i.e.,
+    // [A2 | B2 | c2].
+    // Thus, the coefficients of the parameters after substitution become
+    // (A2 • X + B2)
+    // and the constant terms become
+    // (A2 • y + c2).
+    // The region is therefore given by
+    // (A2 • X + B2) p + (A2 • y + c2) ≥ 0
+    // This is equivalent to A2 • [X | y] + [B2 | c2]
+    // Thus we premultiply [X | y] with each row of A2
+    // and add each row of [B2 | c2].
+    for (unsigned i = 0; i < numIneqs - numVars; i++) {
       activeRegion.setRow(i, (*vertex).preMultiplyWithRow(a2.getRow(i)));
-      activeRegion.addToRow(i, b2c2.getRow(i), Fraction(1, 1));
+      activeRegion.addToRow(i, b2c2.getRow(i), 1);
     }
 
-    activeRegionNorm = IntMatrix(n - d, p + 1);
-    activeRegionRel =
-        IntegerRelation(PresburgerSpace::getRelationSpace(0, p, 0, 0));
-    lcmDenoms = 1;
-    for (unsigned i = 0; i < n - d; i++) {
-      for (unsigned j = 0; j < p + 1; j++)
+    // We convert the representation of the active region to an integers-only
+    // form so as to store it as an PresburgerRelation.
+    // We do this by taking the LCM of the denominators of all the coefficients
+    // and multiplying by it throughout.
+    IntMatrix activeRegionNorm = IntMatrix(numIneqs - numVars, numParams + 1);
+    IntegerRelation activeRegionRel =
+        IntegerRelation(PresburgerSpace::getRelationSpace(0, numParams, 0, 0));
+    MPInt lcmDenoms = MPInt(1);
+    for (unsigned i = 0; i < numIneqs - numVars; i++) {
+      for (unsigned j = 0; j < numParams + 1; j++)
         lcmDenoms = lcm(lcmDenoms, activeRegion(i, j).den);
-      for (unsigned j = 0; j < p + 1; j++)
+      for (unsigned j = 0; j < numParams + 1; j++)
         activeRegionNorm(i, j) =
-            (activeRegion(i, j) * Fraction(lcmDenoms, 1)).getAsInteger();
+            (activeRegion(i, j) * lcmDenoms).getAsInteger();
 
       activeRegionRel.addInequality(activeRegionNorm.getRow(i));
     }
@@ -308,17 +346,45 @@ mlir::presburger::detail::polytopeGeneratingFunction(PolyhedronH poly) {
     activeRegions.push_back(PresburgerRelation(activeRegionRel));
   }
 
-  // Clauss-Loechner chamber decomposition
+  // Now, we use Clauss-Loechner decomposition to identify regions in parameter
+  // space where each vertex is active. These regions (chambers) have the
+  // property that no two of them have a full-dimensional intersection, i.e.,
+  // they may share "faces" or "edges", but their intersection can only have
+  // up to numVars-1 dimensions.
+
+  // We maintain a list of regions and their associated vertex sets,
+  // initialized with the first vertex and its corresponding activity region.
   std::vector<std::pair<PresburgerRelation, std::vector<unsigned>>> chambers = {
       std::make_pair(activeRegions[0], std::vector({0u}))};
+  // Note that instead of storing lists of actual vertices, we store lists
+  // of indices. Thus the set {2, 3, 4} represents the vertex set
+  // {vertices[2], vertices[3], vertices[4]}.
+
   std::vector<std::pair<PresburgerRelation, std::vector<unsigned>>> newChambers;
+
+  // We iterate over the vertex set.
+  // For each vertex v_j and its activity region R_j,
+  // we examine all the current chambers R_i.
+  // If R_j has a full-dimensional intersection with an existing chamber R_i,
+  // then that chamber is replaced by two new ones:
+  // 1. the intersection R_i \cap R_j, where v_j is active;
+  // 2. the difference R_i - R_j, where v_j is inactive.
+  // Once we have examined all R_i, we add a final chamber
+  // R_j - (union of all existing chambers),
+  // in which only v_j is active.
   for (unsigned j = 1; j < vertices.size(); j++) {
     newChambers.clear();
+
     PresburgerRelation r_j = activeRegions[j];
     ParamPoint v_j = vertices[j];
+
     for (unsigned i = 0; i < chambers.size(); i++) {
       auto [r_i, v_i] = chambers[i];
 
+      // First, we check if the intersection of R_j and R_i.
+      // It is a disjoint union of convex regions in the parameter space,
+      // and so we know that it is full-dimensional if any of the disjuncts
+      // is full-dimensional.
       PresburgerRelation intersection = r_i.intersect(r_j);
       bool isFullDim = false;
       for (auto disjunct : intersection.getAllDisjuncts())
@@ -326,10 +392,14 @@ mlir::presburger::detail::polytopeGeneratingFunction(PolyhedronH poly) {
           isFullDim = true;
           break;
         }
-      isFullDim = (p == 0) || isFullDim;
+      isFullDim = (numParams == 0) || isFullDim;
+
+      // If the intersection is not full-dimensional, we do not modify
+      // the chamber list.
       if (!isFullDim)
         newChambers.push_back(chambers[i]);
       else {
+        // If it is, we add the intersection and the difference as new chambers.
         PresburgerRelation subtraction = r_i.subtract(r_j);
         newChambers.push_back(std::make_pair(subtraction, v_i));
 
@@ -337,11 +407,14 @@ mlir::presburger::detail::polytopeGeneratingFunction(PolyhedronH poly) {
         newChambers.push_back(std::make_pair(intersection, v_i));
       }
     }
+
+    // Finally we compute the chamber where only v_j is active by subtracting
+    // all existing chambers from R_j.
     for (auto chamber : newChambers)
       r_j = r_j.subtract(chamber.first);
-
     newChambers.push_back(std::make_pair(r_j, std::vector({j})));
 
+    // We filter `chambers` to remove empty regions.
     chambers.clear();
     for (auto chamber : newChambers) {
       bool empty = true;
@@ -355,18 +428,27 @@ mlir::presburger::detail::polytopeGeneratingFunction(PolyhedronH poly) {
     }
   }
 
-  SmallVector<MPInt> ineq(d + 1);
+  // Now, we compute the generating function. For each chamber, we iterate over
+  // the vertices active in it, and compute the generating function for each
+  // of them. The sum of these generating functions is the GF corresponding to
+  // the entire polytope.
+  SmallVector<MPInt> ineq(numVars + 1);
   for (auto chamber : chambers) {
-    chamberGf = GeneratingFunction(p, {}, {}, {});
+    GeneratingFunction chamberGf(numParams, {}, {}, {});
     for (unsigned i : chamber.second) {
-      tgtCone = defineHRep(d);
-      for (unsigned j = 0; j < d; j++) {
-        for (unsigned k = 0; k < d; k++)
+      // We collect the inequalities corresponding to each vertex.
+      // We only need the coefficients of the variables (NOT the parameters)
+      // as the generating function only depends on these.
+      ConeH tgtCone = defineHRep(numVars);
+      for (unsigned j = 0; j < numVars; j++) {
+        for (unsigned k = 0; k < numVars; k++)
           ineq[k] = subsets[i](j, k);
-        ineq[d] = subsets[i](j, d + p);
+        ineq[numVars] = subsets[i](j, numVars + numParams);
         tgtCone.addInequality(ineq);
       }
-      unimodCones = {std::make_pair(1, tgtCone)};
+      // We assume that the tangent cone is unimodular.
+      SmallVector<std::pair<int, ConeH>, 4> unimodCones = {
+          std::make_pair(1, tgtCone)};
       for (auto signedCone : unimodCones)
         chamberGf =
             chamberGf + unimodularConeGeneratingFunction(

>From 6176e0cc4008946f1ada6bea336d36dbee850ee5 Mon Sep 17 00:00:00 2001
From: Abhinav271828 <abhinav.m at research.iiit.ac.in>
Date: Mon, 22 Jan 2024 20:24:54 +0530
Subject: [PATCH 06/26] Doc

---
 .../mlir/Analysis/Presburger/IntegerRelation.h   |  3 +++
 mlir/lib/Analysis/Presburger/IntegerRelation.cpp | 16 ++++++++++++----
 2 files changed, 15 insertions(+), 4 deletions(-)

diff --git a/mlir/include/mlir/Analysis/Presburger/IntegerRelation.h b/mlir/include/mlir/Analysis/Presburger/IntegerRelation.h
index e7b03c44d33e38c..141dc1f8fa97ed2 100644
--- a/mlir/include/mlir/Analysis/Presburger/IntegerRelation.h
+++ b/mlir/include/mlir/Analysis/Presburger/IntegerRelation.h
@@ -711,8 +711,11 @@ class IntegerRelation {
   /// return `this \ set`.
   PresburgerRelation subtract(const PresburgerRelation &set) const;
 
+  // Remove equalities which have only zero coefficients.
   void removeTrivialEqualities();
 
+  // Verify whether the relation is full-dimensional, i.e.,
+  // has the same number of dimensions as the number of variables.
   bool isFullDim();
 
   void print(raw_ostream &os) const;
diff --git a/mlir/lib/Analysis/Presburger/IntegerRelation.cpp b/mlir/lib/Analysis/Presburger/IntegerRelation.cpp
index edb95b4072e1d85..b05075b2892ace0 100644
--- a/mlir/lib/Analysis/Presburger/IntegerRelation.cpp
+++ b/mlir/lib/Analysis/Presburger/IntegerRelation.cpp
@@ -2511,24 +2511,32 @@ void IntegerRelation::removeTrivialEqualities() {
 }
 
 bool IntegerRelation::isFullDim() {
-  removeTrivialEqualities(); // to implement: remove equalities that are `0 = 0`
+  removeTrivialEqualities();
 
+  // If there is a non-trivial equality, the space cannot be full-dimensional.
   if (getNumEqualities() > 0)
-    return true;
+    return false;
 
+  // If along the direction of any of the inequalities, the upper and lower
+  // optima are the same, then the region is not full-dimensional.
   Simplex simplex(*this);
   for (unsigned i = 0; i < getNumInequalities(); i++) {
     auto ineq = inequalities.getRow(i);
     auto upOpt = simplex.computeOptimum(Simplex::Direction::Up, ineq);
     auto downOpt = simplex.computeOptimum(Simplex::Direction::Down, ineq);
+
     if (upOpt.getKind() == OptimumKind::Unbounded ||
         downOpt.getKind() == OptimumKind::Unbounded)
       continue;
+
+    // Check if the upper and lower optima are equal.
     if (upOpt.getKind() == OptimumKind::Bounded &&
-        downOpt.getKind() == OptimumKind::Bounded &&
-        *upOpt == *downOpt) // ineq == 0 holds for this
+        downOpt.getKind() == OptimumKind::Bounded && *upOpt == *downOpt)
       return false;
   }
+  // If none of the inequalities were such that the upper and lower optima
+  // along their direction were equal, then we conclude that the region is full
+  // dimensional.
   return true;
 }
 

>From 2442f76e7b1a47452cea9f20067a94d0e14df1ae Mon Sep 17 00:00:00 2001
From: Abhinav271828 <abhinav.m at research.iiit.ac.in>
Date: Mon, 22 Jan 2024 20:47:09 +0530
Subject: [PATCH 07/26] ../mlir/lib/Analysis/Presburger/Barvinok.cpp

---
 mlir/lib/Analysis/Presburger/IntegerRelation.cpp | 10 ++++++----
 1 file changed, 6 insertions(+), 4 deletions(-)

diff --git a/mlir/lib/Analysis/Presburger/IntegerRelation.cpp b/mlir/lib/Analysis/Presburger/IntegerRelation.cpp
index b05075b2892ace0..1683f5cf41aee2e 100644
--- a/mlir/lib/Analysis/Presburger/IntegerRelation.cpp
+++ b/mlir/lib/Analysis/Presburger/IntegerRelation.cpp
@@ -2500,9 +2500,9 @@ void IntegerRelation::printSpace(raw_ostream &os) const {
 
 void IntegerRelation::removeTrivialEqualities() {
   bool flag;
-  for (unsigned i = 0; i < getNumEqualities(); i++) {
+  for (unsigned i = 0, e = getNumInequalities(); i < e; i++) {
     flag = true;
-    for (unsigned j = 0; j < getNumVars() + 1; j++)
+    for (unsigned j = 0, f = getNumVars(); j < f + 1; j++)
       if (atEq(i, j) != 0)
         flag = false;
     if (flag)
@@ -2513,14 +2513,16 @@ void IntegerRelation::removeTrivialEqualities() {
 bool IntegerRelation::isFullDim() {
   removeTrivialEqualities();
 
+  unsigned e = getNumInequalities();
+
   // If there is a non-trivial equality, the space cannot be full-dimensional.
-  if (getNumEqualities() > 0)
+  if (e > 0)
     return false;
 
   // If along the direction of any of the inequalities, the upper and lower
   // optima are the same, then the region is not full-dimensional.
   Simplex simplex(*this);
-  for (unsigned i = 0; i < getNumInequalities(); i++) {
+  for (unsigned i = 0; i < e; i++) {
     auto ineq = inequalities.getRow(i);
     auto upOpt = simplex.computeOptimum(Simplex::Direction::Up, ineq);
     auto downOpt = simplex.computeOptimum(Simplex::Direction::Down, ineq);

>From 57ce7b013dce3c49f5fb20464d0ddac35273cd2c Mon Sep 17 00:00:00 2001
From: Abhinav271828 <abhinav.m at research.iiit.ac.in>
Date: Mon, 22 Jan 2024 20:47:18 +0530
Subject: [PATCH 08/26] Simplifications

---
 mlir/lib/Analysis/Presburger/Barvinok.cpp | 49 ++++++++++-------------
 1 file changed, 22 insertions(+), 27 deletions(-)

diff --git a/mlir/lib/Analysis/Presburger/Barvinok.cpp b/mlir/lib/Analysis/Presburger/Barvinok.cpp
index 3afc9ad8638ccd2..1a3dd9a32e50af8 100644
--- a/mlir/lib/Analysis/Presburger/Barvinok.cpp
+++ b/mlir/lib/Analysis/Presburger/Barvinok.cpp
@@ -263,12 +263,14 @@ mlir::presburger::detail::polytopeGeneratingFunction(PolyhedronH poly) {
   FracMatrix b2c2(numIneqs - numVars, numParams + 1);
 
   // We iterate over all subsets of inequalities with cardinality numVars,
-  // using bitsets up to 2^numIneqs to enumerate.
+  // using bitsets to enumerate.
+  // The largest possible bitset that corresponds to such a subset can be
+  // written as numVar 1's followed by (numIneqs - numVars) 0's.
+  unsigned upperBound = ((1ul << numVars) - 1ul)
+           << (numIneqs - numVars);
   for (std::bitset<16> indicator(((1ul << numVars) - 1ul)
                                  << (numIneqs - numVars));
-       indicator.to_ulong() <=
-       ((1ul << numVars) - 1ul)
-           << (numIneqs - numVars); // d 1's followed by n-numVars 0's
+       indicator.to_ulong() <= upperBound;
        indicator = std::bitset<16>(indicator.to_ulong() - 1)) {
 
     if (indicator.count() != numVars)
@@ -372,13 +374,13 @@ mlir::presburger::detail::polytopeGeneratingFunction(PolyhedronH poly) {
   // Once we have examined all R_i, we add a final chamber
   // R_j - (union of all existing chambers),
   // in which only v_j is active.
-  for (unsigned j = 1; j < vertices.size(); j++) {
+  for (unsigned j = 1, e = vertices.size(); j < e; j++) {
     newChambers.clear();
 
     PresburgerRelation r_j = activeRegions[j];
     ParamPoint v_j = vertices[j];
 
-    for (unsigned i = 0; i < chambers.size(); i++) {
+    for (unsigned i = 0, f = chambers.size(); i < f; i++) {
       auto [r_i, v_i] = chambers[i];
 
       // First, we check if the intersection of R_j and R_i.
@@ -386,13 +388,7 @@ mlir::presburger::detail::polytopeGeneratingFunction(PolyhedronH poly) {
       // and so we know that it is full-dimensional if any of the disjuncts
       // is full-dimensional.
       PresburgerRelation intersection = r_i.intersect(r_j);
-      bool isFullDim = false;
-      for (auto disjunct : intersection.getAllDisjuncts())
-        if (disjunct.isFullDim()) {
-          isFullDim = true;
-          break;
-        }
-      isFullDim = (numParams == 0) || isFullDim;
+      bool isFullDim = numParams == 0 || llvm::any_of(intersection.getAllDisjuncts(), [&](IntegerRelation disjunct) -> bool { return disjunct.isFullDim(); });
 
       // If the intersection is not full-dimensional, we do not modify
       // the chamber list.
@@ -410,20 +406,16 @@ mlir::presburger::detail::polytopeGeneratingFunction(PolyhedronH poly) {
 
     // Finally we compute the chamber where only v_j is active by subtracting
     // all existing chambers from R_j.
-    for (auto chamber : newChambers)
+    for (const std::pair<PresburgerRelation, std::vector<unsigned>> &chamber : newChambers)
       r_j = r_j.subtract(chamber.first);
     newChambers.push_back(std::make_pair(r_j, std::vector({j})));
 
     // We filter `chambers` to remove empty regions.
     chambers.clear();
-    for (auto chamber : newChambers) {
-      bool empty = true;
-      for (auto disjunct : chamber.first.getAllDisjuncts())
-        if (!disjunct.isEmpty()) {
-          empty = false;
-          break;
-        }
-      if (!empty)
+    for (const std::pair<PresburgerRelation, std::vector<unsigned>> &chamber : newChambers) {
+      auto [r, v] = chamber;
+      bool isEmpty = llvm::all_of(r.getAllDisjuncts(), [&](IntegerRelation disjunct) -> bool { return disjunct.isEmpty(); });
+      if (!isEmpty)
         chambers.push_back(chamber);
     }
   }
@@ -433,9 +425,10 @@ mlir::presburger::detail::polytopeGeneratingFunction(PolyhedronH poly) {
   // of them. The sum of these generating functions is the GF corresponding to
   // the entire polytope.
   SmallVector<MPInt> ineq(numVars + 1);
-  for (auto chamber : chambers) {
+  for (const std::pair<PresburgerRelation, std::vector<unsigned>> &chamber : chambers) {
+    auto [region_j, vertices_j] = chamber;
     GeneratingFunction chamberGf(numParams, {}, {}, {});
-    for (unsigned i : chamber.second) {
+    for (unsigned i : vertices_j) {
       // We collect the inequalities corresponding to each vertex.
       // We only need the coefficients of the variables (NOT the parameters)
       // as the generating function only depends on these.
@@ -449,12 +442,14 @@ mlir::presburger::detail::polytopeGeneratingFunction(PolyhedronH poly) {
       // We assume that the tangent cone is unimodular.
       SmallVector<std::pair<int, ConeH>, 4> unimodCones = {
           std::make_pair(1, tgtCone)};
-      for (auto signedCone : unimodCones)
+      for (std::pair<int, ConeH> signedCone : unimodCones) {
+        auto [sign, cone] = signedCone;
         chamberGf =
             chamberGf + unimodularConeGeneratingFunction(
-                            vertices[i], signedCone.first, signedCone.second);
+                            vertices[i], sign, cone);
+      }
     }
-    gf.push_back(std::make_pair(chamber.first, chamberGf));
+    gf.push_back(std::make_pair(region_j, chamberGf));
   }
   return gf;
 }

>From 09a89d7018ee8279d4abf5724636a5180dd07c09 Mon Sep 17 00:00:00 2001
From: Abhinav271828 <abhinav.m at research.iiit.ac.in>
Date: Mon, 22 Jan 2024 20:50:00 +0530
Subject: [PATCH 09/26] Formatting

---
 mlir/lib/Analysis/Presburger/Barvinok.cpp | 27 ++++++++++++++---------
 1 file changed, 17 insertions(+), 10 deletions(-)

diff --git a/mlir/lib/Analysis/Presburger/Barvinok.cpp b/mlir/lib/Analysis/Presburger/Barvinok.cpp
index 1a3dd9a32e50af8..b3e68d1f65f2561 100644
--- a/mlir/lib/Analysis/Presburger/Barvinok.cpp
+++ b/mlir/lib/Analysis/Presburger/Barvinok.cpp
@@ -266,8 +266,7 @@ mlir::presburger::detail::polytopeGeneratingFunction(PolyhedronH poly) {
   // using bitsets to enumerate.
   // The largest possible bitset that corresponds to such a subset can be
   // written as numVar 1's followed by (numIneqs - numVars) 0's.
-  unsigned upperBound = ((1ul << numVars) - 1ul)
-           << (numIneqs - numVars);
+  unsigned upperBound = ((1ul << numVars) - 1ul) << (numIneqs - numVars);
   for (std::bitset<16> indicator(((1ul << numVars) - 1ul)
                                  << (numIneqs - numVars));
        indicator.to_ulong() <= upperBound;
@@ -388,7 +387,11 @@ mlir::presburger::detail::polytopeGeneratingFunction(PolyhedronH poly) {
       // and so we know that it is full-dimensional if any of the disjuncts
       // is full-dimensional.
       PresburgerRelation intersection = r_i.intersect(r_j);
-      bool isFullDim = numParams == 0 || llvm::any_of(intersection.getAllDisjuncts(), [&](IntegerRelation disjunct) -> bool { return disjunct.isFullDim(); });
+      bool isFullDim =
+          numParams == 0 || llvm::any_of(intersection.getAllDisjuncts(),
+                                         [&](IntegerRelation disjunct) -> bool {
+                                           return disjunct.isFullDim();
+                                         });
 
       // If the intersection is not full-dimensional, we do not modify
       // the chamber list.
@@ -406,15 +409,19 @@ mlir::presburger::detail::polytopeGeneratingFunction(PolyhedronH poly) {
 
     // Finally we compute the chamber where only v_j is active by subtracting
     // all existing chambers from R_j.
-    for (const std::pair<PresburgerRelation, std::vector<unsigned>> &chamber : newChambers)
+    for (const std::pair<PresburgerRelation, std::vector<unsigned>> &chamber :
+         newChambers)
       r_j = r_j.subtract(chamber.first);
     newChambers.push_back(std::make_pair(r_j, std::vector({j})));
 
     // We filter `chambers` to remove empty regions.
     chambers.clear();
-    for (const std::pair<PresburgerRelation, std::vector<unsigned>> &chamber : newChambers) {
+    for (const std::pair<PresburgerRelation, std::vector<unsigned>> &chamber :
+         newChambers) {
       auto [r, v] = chamber;
-      bool isEmpty = llvm::all_of(r.getAllDisjuncts(), [&](IntegerRelation disjunct) -> bool { return disjunct.isEmpty(); });
+      bool isEmpty = llvm::all_of(
+          r.getAllDisjuncts(),
+          [&](IntegerRelation disjunct) -> bool { return disjunct.isEmpty(); });
       if (!isEmpty)
         chambers.push_back(chamber);
     }
@@ -425,7 +432,8 @@ mlir::presburger::detail::polytopeGeneratingFunction(PolyhedronH poly) {
   // of them. The sum of these generating functions is the GF corresponding to
   // the entire polytope.
   SmallVector<MPInt> ineq(numVars + 1);
-  for (const std::pair<PresburgerRelation, std::vector<unsigned>> &chamber : chambers) {
+  for (const std::pair<PresburgerRelation, std::vector<unsigned>> &chamber :
+       chambers) {
     auto [region_j, vertices_j] = chamber;
     GeneratingFunction chamberGf(numParams, {}, {}, {});
     for (unsigned i : vertices_j) {
@@ -444,9 +452,8 @@ mlir::presburger::detail::polytopeGeneratingFunction(PolyhedronH poly) {
           std::make_pair(1, tgtCone)};
       for (std::pair<int, ConeH> signedCone : unimodCones) {
         auto [sign, cone] = signedCone;
-        chamberGf =
-            chamberGf + unimodularConeGeneratingFunction(
-                            vertices[i], sign, cone);
+        chamberGf = chamberGf +
+                    unimodularConeGeneratingFunction(vertices[i], sign, cone);
       }
     }
     gf.push_back(std::make_pair(region_j, chamberGf));

>From 2562fefb9a4290182b96a3ca654da5cad83f55d1 Mon Sep 17 00:00:00 2001
From: Abhinav271828 <abhinav.m at research.iiit.ac.in>
Date: Tue, 23 Jan 2024 02:43:33 +0530
Subject: [PATCH 10/26] Doc

---
 .../mlir/Analysis/Presburger/Barvinok.h       |   4 +-
 .../Analysis/Presburger/IntegerRelation.h     |   2 +-
 .../mlir/Analysis/Presburger/Simplex.h        |   4 +
 mlir/lib/Analysis/Presburger/Barvinok.cpp     | 106 ++++++++++--------
 mlir/lib/Analysis/Presburger/Simplex.cpp      |   8 ++
 5 files changed, 74 insertions(+), 50 deletions(-)

diff --git a/mlir/include/mlir/Analysis/Presburger/Barvinok.h b/mlir/include/mlir/Analysis/Presburger/Barvinok.h
index a5bc49d79408c6c..b4d1d6de5307978 100644
--- a/mlir/include/mlir/Analysis/Presburger/Barvinok.h
+++ b/mlir/include/mlir/Analysis/Presburger/Barvinok.h
@@ -88,8 +88,8 @@ GeneratingFunction unimodularConeGeneratingFunction(ParamPoint vertex, int sign,
 /// Find the solution of a set of equations that express affine constraints
 /// between a set of variables and a set of parameters. The solution expresses
 /// each variable as an affine function of the parameters.
-/// If there is no solution, return std::null.
-std::optional<ParamPoint> findVertex(IntMatrix equations);
+/// If there is no solution, return null.
+std::optional<ParamPoint> solveParametricEquations(FracMatrix equations);
 
 /// Compute the generating function corresponding to a polytope.
 /// All tangent cones of the polytope must be unimodular.
diff --git a/mlir/include/mlir/Analysis/Presburger/IntegerRelation.h b/mlir/include/mlir/Analysis/Presburger/IntegerRelation.h
index 141dc1f8fa97ed2..b60d68aab4cf405 100644
--- a/mlir/include/mlir/Analysis/Presburger/IntegerRelation.h
+++ b/mlir/include/mlir/Analysis/Presburger/IntegerRelation.h
@@ -715,7 +715,7 @@ class IntegerRelation {
   void removeTrivialEqualities();
 
   // Verify whether the relation is full-dimensional, i.e.,
-  // has the same number of dimensions as the number of variables.
+  // no equality holds for the relation.
   bool isFullDim();
 
   void print(raw_ostream &os) const;
diff --git a/mlir/include/mlir/Analysis/Presburger/Simplex.h b/mlir/include/mlir/Analysis/Presburger/Simplex.h
index 9482f69b31cd666..1939a41e52bc6d1 100644
--- a/mlir/include/mlir/Analysis/Presburger/Simplex.h
+++ b/mlir/include/mlir/Analysis/Presburger/Simplex.h
@@ -771,6 +771,10 @@ class Simplex : public SimplexBase {
   std::pair<MaybeOptimum<MPInt>, MaybeOptimum<MPInt>>
   computeIntegerBounds(ArrayRef<MPInt> coeffs);
 
+  /// Returns false if the given equality has distinct upper and lower bounds
+  /// in the simplex.
+  bool isValidEquality(const ArrayRef<MPInt> coeffs);
+
   /// Returns true if the polytope is unbounded, i.e., extends to infinity in
   /// some direction. Otherwise, returns false.
   bool isUnbounded();
diff --git a/mlir/lib/Analysis/Presburger/Barvinok.cpp b/mlir/lib/Analysis/Presburger/Barvinok.cpp
index b3e68d1f65f2561..6452b2d37838dc1 100644
--- a/mlir/lib/Analysis/Presburger/Barvinok.cpp
+++ b/mlir/lib/Analysis/Presburger/Barvinok.cpp
@@ -158,65 +158,70 @@ GeneratingFunction mlir::presburger::detail::unimodularConeGeneratingFunction(
 /// therefore represented as a matrix of size d x (p+1).
 /// If there is no solution, we return null.
 std::optional<ParamPoint>
-mlir::presburger::detail::findVertex(IntMatrix equations) {
+mlir::presburger::detail::solveParametricEquations(FracMatrix equations) {
   // equations is a d x (d + p + 1) matrix.
   // Each row represents an equation.
-
   unsigned numEqs = equations.getNumRows();
   unsigned numCols = equations.getNumColumns();
 
   // First, we check that the system has a solution, and return
   // null if not.
-  IntMatrix coeffs(numEqs, numEqs);
+  FracMatrix coeffs(numEqs, numEqs);
   for (unsigned i = 0; i < numEqs; i++)
     for (unsigned j = 0; j < numEqs; j++)
       coeffs(i, j) = equations(i, j);
 
+  // If the determinant is zero, there is no unique solution.
+  // Thus we return null.
   if (coeffs.determinant() == 0)
     return std::nullopt;
 
-  // We work with rational numbers.
-  FracMatrix equationsF(equations);
-
   for (unsigned i = 0; i < numEqs; ++i) {
     // First ensure that the diagonal element is nonzero, by swapping
-    // it with a nonzero row.
-    if (equationsF(i, i) == 0) {
+    // it with a row that is non-zero at column i.
+    if (equations(i, i) == 0) {
       for (unsigned j = i + 1; j < numEqs; ++j) {
-        if (equationsF(j, i) != 0) {
-          equationsF.swapRows(j, i);
-          break;
-        }
+        if (equations(j, i) == 0)
+          continue;
+        equations.swapRows(j, i);
+        break;
       }
     }
 
-    Fraction b = equationsF(i, i);
+    Fraction diagElement = equations(i, i);
 
-    // Set all elements except the diagonal to zero.
+    // Apply row operations to make all elements except the diagonal to zero.
     for (unsigned j = 0; j < numEqs; ++j) {
-      if (equationsF(j, i) == 0 || j == i)
+      if (i == j)
         continue;
-      // Set element (j, i) to zero
-      // by subtracting the ith row,
-      // appropriately scaled.
-      Fraction a = equationsF(j, i);
-      equationsF.addToRow(j, equationsF.getRow(i), -a / b);
+      if (equations(j, i) == 0)
+        continue;
+      // Apply row operations to make element (j, i) zero by subtracting the
+      // ith row, appropriately scaled.
+      Fraction currentElement = equations(j, i);
+      equations.addToRow(j, equations.getRow(i), -currentElement / diagElement);
     }
   }
 
   // Rescale diagonal elements to 1.
   for (unsigned i = 0; i < numEqs; ++i) {
-    Fraction a = equationsF(i, i);
+    Fraction diagElement = equations(i, i);
     for (unsigned j = 0; j < numCols; ++j)
-      equationsF(i, j) = equationsF(i, j) / a;
+      equations(i, j) = equations(i, j) / diagElement;
   }
 
-  // We copy the last p+1 columns of the matrix as the values of x_i.
-  // We shift the parameter terms to the RHS, and so flip their sign.
+  // Now we have reduced the equations to the form
+  // x_i + b_1' m_1 + ... + b_p' m_p + c' = 0
+  // i.e. each variable appears exactly once in the system, and has coefficient
+  // one.
+  // Thus we have
+  // x_i = - b_1' m_1 - ... - b_p' m_p - c
+  // and so we return the negation of the last p + 1 columns of the matrix.
+  // We copy these columns and return them.
   ParamPoint vertex(numEqs, numCols - numEqs);
   for (unsigned i = 0; i < numEqs; ++i)
     for (unsigned j = 0; j < numCols - numEqs; ++j)
-      vertex(i, j) = -equationsF(i, numEqs + j);
+      vertex(i, j) = -equations(i, numEqs + j);
 
   return vertex;
 }
@@ -225,7 +230,7 @@ mlir::presburger::detail::findVertex(IntMatrix equations) {
 /// function corresponding to the number of lattice points present. This
 /// algorithm has three main steps:
 /// 1. Enumerate the vertices, by iterating over subsets of inequalities and
-///    checking for solubility.
+///    checking for satisfiability.
 /// 2. For each vertex, identify the tangent cone and compute the generating
 ///    function corresponding to it. The sum of these GFs is the GF of the
 ///    polytope.
@@ -239,17 +244,17 @@ mlir::presburger::detail::findVertex(IntMatrix equations) {
 std::vector<std::pair<PresburgerRelation, GeneratingFunction>>
 mlir::presburger::detail::polytopeGeneratingFunction(PolyhedronH poly) {
   unsigned numVars = poly.getNumRangeVars();
-  unsigned numParams = poly.getNumSymbolVars();
+  unsigned numSymbols = poly.getNumSymbolVars();
   unsigned numIneqs = poly.getNumInequalities();
 
   // The generating function of the polytope is computed as a set of generating
   // functions, each one associated with a region in parameter space (chamber).
   std::vector<std::pair<PresburgerRelation, GeneratingFunction>> gf({});
 
-  // The active region will be defined as activeRegionCoeffs @ p +
+  // The active region will be defined as activeRegionCoeffs x p +
   // activeRegionConstant ≥ 0. The active region is a polyhedron in parameter
   // space.
-  FracMatrix activeRegion(numIneqs - numVars, numParams + 1);
+  FracMatrix activeRegion(numIneqs - numVars, numSymbols + 1);
 
   // These vectors store lists of
   // subsets of inequalities,
@@ -260,12 +265,17 @@ mlir::presburger::detail::polytopeGeneratingFunction(PolyhedronH poly) {
   std::vector<PresburgerRelation> activeRegions;
 
   FracMatrix a2(numIneqs - numVars, numVars);
-  FracMatrix b2c2(numIneqs - numVars, numParams + 1);
+  FracMatrix b2c2(numIneqs - numVars, numSymbols + 1);
 
   // We iterate over all subsets of inequalities with cardinality numVars,
-  // using bitsets to enumerate.
+  // using bitsets of numIneqs bits to enumerate.
+  // For a given set of numIneqs bits, we consider a subset which contains
+  // the i'th inequality if the i'th bit in the bitset is set.
   // The largest possible bitset that corresponds to such a subset can be
   // written as numVar 1's followed by (numIneqs - numVars) 0's.
+  // We start with this and count downwards (in binary), considering only
+  // those numbers whose binary representation has numVars 1's and the
+  // rest 0's.
   unsigned upperBound = ((1ul << numVars) - 1ul) << (numIneqs - numVars);
   for (std::bitset<16> indicator(((1ul << numVars) - 1ul)
                                  << (numIneqs - numVars));
@@ -276,9 +286,9 @@ mlir::presburger::detail::polytopeGeneratingFunction(PolyhedronH poly) {
       continue;
 
     // Collect the inequalities corresponding to the bits which are set.
-    IntMatrix subset(numVars, numVars + numParams + 1);
+    IntMatrix subset(numVars, numVars + numSymbols + 1);
     unsigned j1 = 0, j2 = 0;
-    for (unsigned i = 0; i < numIneqs; i++)
+    for (unsigned i = 0; i < numIneqs; i++) {
       if (indicator.test(i))
         subset.setRow(j1++, poly.getInequality(i));
 
@@ -289,16 +299,18 @@ mlir::presburger::detail::polytopeGeneratingFunction(PolyhedronH poly) {
         // b2c2 stores the coefficients of the parameters and the constant term.
         for (unsigned k = 0; k < numVars; k++)
           a2(j2, k) = poly.atIneq(i, k);
-        for (unsigned k = numVars; k < numVars + numParams + 1; k++)
+        for (unsigned k = numVars; k < numVars + numSymbols + 1; k++)
           b2c2(j2, k - numVars) = poly.atIneq(i, k);
         j2++;
       }
+    }
 
     // Find the vertex, if any, corresponding to the current subset of
     // inequalities.
-    std::optional<ParamPoint> vertex = findVertex(subset); // d x (p+1)
+    std::optional<ParamPoint> vertex =
+        solveParametricEquations(FracMatrix(subset)); // d x (p+1)
 
-    if (vertex == std::nullopt)
+    if (!vertex)
       continue;
     // If this subset corresponds to a vertex, store it.
     vertices.push_back(*vertex);
@@ -330,14 +342,14 @@ mlir::presburger::detail::polytopeGeneratingFunction(PolyhedronH poly) {
     // form so as to store it as an PresburgerRelation.
     // We do this by taking the LCM of the denominators of all the coefficients
     // and multiplying by it throughout.
-    IntMatrix activeRegionNorm = IntMatrix(numIneqs - numVars, numParams + 1);
+    IntMatrix activeRegionNorm = IntMatrix(numIneqs - numVars, numSymbols + 1);
     IntegerRelation activeRegionRel =
-        IntegerRelation(PresburgerSpace::getRelationSpace(0, numParams, 0, 0));
+        IntegerRelation(PresburgerSpace::getRelationSpace(0, numSymbols, 0, 0));
     MPInt lcmDenoms = MPInt(1);
     for (unsigned i = 0; i < numIneqs - numVars; i++) {
-      for (unsigned j = 0; j < numParams + 1; j++)
+      for (unsigned j = 0; j < numSymbols + 1; j++)
         lcmDenoms = lcm(lcmDenoms, activeRegion(i, j).den);
-      for (unsigned j = 0; j < numParams + 1; j++)
+      for (unsigned j = 0; j < numSymbols + 1; j++)
         activeRegionNorm(i, j) =
             (activeRegion(i, j) * lcmDenoms).getAsInteger();
 
@@ -387,11 +399,11 @@ mlir::presburger::detail::polytopeGeneratingFunction(PolyhedronH poly) {
       // and so we know that it is full-dimensional if any of the disjuncts
       // is full-dimensional.
       PresburgerRelation intersection = r_i.intersect(r_j);
-      bool isFullDim =
-          numParams == 0 || llvm::any_of(intersection.getAllDisjuncts(),
-                                         [&](IntegerRelation disjunct) -> bool {
-                                           return disjunct.isFullDim();
-                                         });
+      bool isFullDim = numSymbols == 0 ||
+                       llvm::any_of(intersection.getAllDisjuncts(),
+                                    [&](IntegerRelation disjunct) -> bool {
+                                      return disjunct.isFullDim();
+                                    });
 
       // If the intersection is not full-dimensional, we do not modify
       // the chamber list.
@@ -435,7 +447,7 @@ mlir::presburger::detail::polytopeGeneratingFunction(PolyhedronH poly) {
   for (const std::pair<PresburgerRelation, std::vector<unsigned>> &chamber :
        chambers) {
     auto [region_j, vertices_j] = chamber;
-    GeneratingFunction chamberGf(numParams, {}, {}, {});
+    GeneratingFunction chamberGf(numSymbols, {}, {}, {});
     for (unsigned i : vertices_j) {
       // We collect the inequalities corresponding to each vertex.
       // We only need the coefficients of the variables (NOT the parameters)
@@ -444,7 +456,7 @@ mlir::presburger::detail::polytopeGeneratingFunction(PolyhedronH poly) {
       for (unsigned j = 0; j < numVars; j++) {
         for (unsigned k = 0; k < numVars; k++)
           ineq[k] = subsets[i](j, k);
-        ineq[numVars] = subsets[i](j, numVars + numParams);
+        ineq[numVars] = subsets[i](j, numVars + numSymbols);
         tgtCone.addInequality(ineq);
       }
       // We assume that the tangent cone is unimodular.
diff --git a/mlir/lib/Analysis/Presburger/Simplex.cpp b/mlir/lib/Analysis/Presburger/Simplex.cpp
index 42bbc3363d5830f..f6d24b3f4dbd40e 100644
--- a/mlir/lib/Analysis/Presburger/Simplex.cpp
+++ b/mlir/lib/Analysis/Presburger/Simplex.cpp
@@ -2104,6 +2104,14 @@ Simplex::computeIntegerBounds(ArrayRef<MPInt> coeffs) {
   return {minRoundedUp, maxRoundedDown};
 }
 
+bool Simplex::isValidEquality(ArrayRef<MPInt> coeffs) {
+  auto [downOpt, upOpt] = Simplex::computeIntegerBounds(coeffs);
+  if (upOpt.getKind() == OptimumKind::Bounded &&
+      downOpt.getKind() == OptimumKind::Bounded && *upOpt == *downOpt)
+    return false;
+  return true;
+}
+
 void SimplexBase::print(raw_ostream &os) const {
   os << "rows = " << getNumRows() << ", columns = " << getNumColumns() << "\n";
   if (empty)

>From 90badd7d6786f0f9f2ad4459c0717eda5dac1d5e Mon Sep 17 00:00:00 2001
From: Abhinav271828 <abhinav.m at research.iiit.ac.in>
Date: Tue, 23 Jan 2024 03:00:44 +0530
Subject: [PATCH 11/26] Implement and use getSubMatrix()

---
 .../include/mlir/Analysis/Presburger/Matrix.h |  6 ++++++
 mlir/lib/Analysis/Presburger/Barvinok.cpp     | 19 ++++++++-----------
 mlir/lib/Analysis/Presburger/Matrix.cpp       | 13 +++++++++++++
 3 files changed, 27 insertions(+), 11 deletions(-)

diff --git a/mlir/include/mlir/Analysis/Presburger/Matrix.h b/mlir/include/mlir/Analysis/Presburger/Matrix.h
index 0d4a593a95b1c96..4a5bba17e204531 100644
--- a/mlir/include/mlir/Analysis/Presburger/Matrix.h
+++ b/mlir/include/mlir/Analysis/Presburger/Matrix.h
@@ -184,6 +184,12 @@ class Matrix {
   // Transpose the matrix without modifying it.
   Matrix<T> transpose() const;
 
+  // Copy the cells in the intersection of
+  // the rows between `fromRows` and `toRows` and
+  // the columns between `fromColumns` and `toColumns`, both inclusive.
+  Matrix<T> getSubMatrix(unsigned fromRow, unsigned toRow, unsigned fromColumn,
+                         unsigned toColumn) const;
+
   /// Print the matrix.
   void print(raw_ostream &os) const;
   void dump() const;
diff --git a/mlir/lib/Analysis/Presburger/Barvinok.cpp b/mlir/lib/Analysis/Presburger/Barvinok.cpp
index 6452b2d37838dc1..67af0eb452bf3f1 100644
--- a/mlir/lib/Analysis/Presburger/Barvinok.cpp
+++ b/mlir/lib/Analysis/Presburger/Barvinok.cpp
@@ -164,16 +164,12 @@ mlir::presburger::detail::solveParametricEquations(FracMatrix equations) {
   unsigned numEqs = equations.getNumRows();
   unsigned numCols = equations.getNumColumns();
 
-  // First, we check that the system has a solution, and return
-  // null if not.
-  FracMatrix coeffs(numEqs, numEqs);
-  for (unsigned i = 0; i < numEqs; i++)
-    for (unsigned j = 0; j < numEqs; j++)
-      coeffs(i, j) = equations(i, j);
-
   // If the determinant is zero, there is no unique solution.
   // Thus we return null.
-  if (coeffs.determinant() == 0)
+  if (FracMatrix(equations.getSubMatrix(/*fromRow=*/0, /*toRow=*/numEqs - 1,
+                                        /*fromColumn=*/0,
+                                        /*toColumn=*/numEqs - 1))
+          .determinant() == 0)
     return std::nullopt;
 
   for (unsigned i = 0; i < numEqs; ++i) {
@@ -218,10 +214,11 @@ mlir::presburger::detail::solveParametricEquations(FracMatrix equations) {
   // x_i = - b_1' m_1 - ... - b_p' m_p - c
   // and so we return the negation of the last p + 1 columns of the matrix.
   // We copy these columns and return them.
-  ParamPoint vertex(numEqs, numCols - numEqs);
+  ParamPoint vertex =
+      equations.getSubMatrix(/*fromRow=*/0, /*toRow=*/numEqs - 1,
+                             /*fromRow=*/numEqs, /*toRow=*/numCols);
   for (unsigned i = 0; i < numEqs; ++i)
-    for (unsigned j = 0; j < numCols - numEqs; ++j)
-      vertex(i, j) = -equations(i, numEqs + j);
+    vertex.negateRow(i);
 
   return vertex;
 }
diff --git a/mlir/lib/Analysis/Presburger/Matrix.cpp b/mlir/lib/Analysis/Presburger/Matrix.cpp
index bd7f7f58a932f3e..be8dac636e0af61 100644
--- a/mlir/lib/Analysis/Presburger/Matrix.cpp
+++ b/mlir/lib/Analysis/Presburger/Matrix.cpp
@@ -354,6 +354,19 @@ static void modEntryColumnOperation(Matrix<MPInt> &m, unsigned row,
   otherMatrix.addToColumn(sourceCol, targetCol, ratio);
 }
 
+template <typename T>
+Matrix<T> Matrix<T>::getSubMatrix(unsigned fromRow, unsigned toRow,
+                                  unsigned fromColumn,
+                                  unsigned toColumn) const {
+  assert(toRow >= fromRow && "end of row range must be after beginning!");
+  assert(toColumn >= fromColumn && "end of row range must be after beginning!");
+  Matrix<T> subMatrix(toRow - fromRow + 1, toColumn - fromColumn + 1);
+  for (unsigned i = fromRow; i <= toRow; i++)
+    for (unsigned j = fromColumn; j <= toColumn; j++)
+      subMatrix(i, j) = at(i, j);
+  return subMatrix;
+}
+
 template <typename T>
 void Matrix<T>::print(raw_ostream &os) const {
   for (unsigned row = 0; row < nRows; ++row) {

>From 7ab464a315e7d85da57a776820f8cd04560d3752 Mon Sep 17 00:00:00 2001
From: Abhinav271828 <abhinav.m at research.iiit.ac.in>
Date: Tue, 23 Jan 2024 03:14:30 +0530
Subject: [PATCH 12/26] Fix names

---
 mlir/lib/Analysis/Presburger/Barvinok.cpp | 26 +++++++++++------------
 1 file changed, 13 insertions(+), 13 deletions(-)

diff --git a/mlir/lib/Analysis/Presburger/Barvinok.cpp b/mlir/lib/Analysis/Presburger/Barvinok.cpp
index 67af0eb452bf3f1..8c4de595b685c49 100644
--- a/mlir/lib/Analysis/Presburger/Barvinok.cpp
+++ b/mlir/lib/Analysis/Presburger/Barvinok.cpp
@@ -385,17 +385,17 @@ mlir::presburger::detail::polytopeGeneratingFunction(PolyhedronH poly) {
   for (unsigned j = 1, e = vertices.size(); j < e; j++) {
     newChambers.clear();
 
-    PresburgerRelation r_j = activeRegions[j];
-    ParamPoint v_j = vertices[j];
+    PresburgerRelation newRegion = activeRegions[j];
+    ParamPoint newVertex = vertices[j];
 
     for (unsigned i = 0, f = chambers.size(); i < f; i++) {
-      auto [r_i, v_i] = chambers[i];
+      auto [currentRegion, currentVertices] = chambers[i];
 
       // First, we check if the intersection of R_j and R_i.
       // It is a disjoint union of convex regions in the parameter space,
       // and so we know that it is full-dimensional if any of the disjuncts
       // is full-dimensional.
-      PresburgerRelation intersection = r_i.intersect(r_j);
+      PresburgerRelation intersection = currentRegion.intersect(newRegion);
       bool isFullDim = numSymbols == 0 ||
                        llvm::any_of(intersection.getAllDisjuncts(),
                                     [&](IntegerRelation disjunct) -> bool {
@@ -408,11 +408,11 @@ mlir::presburger::detail::polytopeGeneratingFunction(PolyhedronH poly) {
         newChambers.push_back(chambers[i]);
       else {
         // If it is, we add the intersection and the difference as new chambers.
-        PresburgerRelation subtraction = r_i.subtract(r_j);
-        newChambers.push_back(std::make_pair(subtraction, v_i));
+        PresburgerRelation subtraction = currentRegion.subtract(newRegion);
+        newChambers.push_back(std::make_pair(subtraction, currentVertices));
 
-        v_i.push_back(j);
-        newChambers.push_back(std::make_pair(intersection, v_i));
+        currentVertices.push_back(j);
+        newChambers.push_back(std::make_pair(intersection, currentVertices));
       }
     }
 
@@ -420,8 +420,8 @@ mlir::presburger::detail::polytopeGeneratingFunction(PolyhedronH poly) {
     // all existing chambers from R_j.
     for (const std::pair<PresburgerRelation, std::vector<unsigned>> &chamber :
          newChambers)
-      r_j = r_j.subtract(chamber.first);
-    newChambers.push_back(std::make_pair(r_j, std::vector({j})));
+      newRegion = newRegion.subtract(chamber.first);
+    newChambers.push_back(std::make_pair(newRegion, std::vector({j})));
 
     // We filter `chambers` to remove empty regions.
     chambers.clear();
@@ -443,9 +443,9 @@ mlir::presburger::detail::polytopeGeneratingFunction(PolyhedronH poly) {
   SmallVector<MPInt> ineq(numVars + 1);
   for (const std::pair<PresburgerRelation, std::vector<unsigned>> &chamber :
        chambers) {
-    auto [region_j, vertices_j] = chamber;
+    auto [currentRegion, currentVertices] = chamber;
     GeneratingFunction chamberGf(numSymbols, {}, {}, {});
-    for (unsigned i : vertices_j) {
+    for (unsigned i : currentVertices) {
       // We collect the inequalities corresponding to each vertex.
       // We only need the coefficients of the variables (NOT the parameters)
       // as the generating function only depends on these.
@@ -465,7 +465,7 @@ mlir::presburger::detail::polytopeGeneratingFunction(PolyhedronH poly) {
                     unimodularConeGeneratingFunction(vertices[i], sign, cone);
       }
     }
-    gf.push_back(std::make_pair(region_j, chamberGf));
+    gf.push_back(std::make_pair(currentRegion, chamberGf));
   }
   return gf;
 }

>From c6afaa219bf8c01bf924887fca7bacdf0ccfbf43 Mon Sep 17 00:00:00 2001
From: Abhinav271828 <abhinav.m at research.iiit.ac.in>
Date: Tue, 23 Jan 2024 17:22:56 +0530
Subject: [PATCH 13/26] Factor out functions

---
 .../mlir/Analysis/Presburger/Barvinok.h       |  14 ++
 .../include/mlir/Analysis/Presburger/Matrix.h |  11 +
 mlir/lib/Analysis/Presburger/Barvinok.cpp     | 206 +++++++++---------
 mlir/lib/Analysis/Presburger/Matrix.cpp       |  32 +++
 4 files changed, 157 insertions(+), 106 deletions(-)

diff --git a/mlir/include/mlir/Analysis/Presburger/Barvinok.h b/mlir/include/mlir/Analysis/Presburger/Barvinok.h
index b4d1d6de5307978..34055e93242c4ca 100644
--- a/mlir/include/mlir/Analysis/Presburger/Barvinok.h
+++ b/mlir/include/mlir/Analysis/Presburger/Barvinok.h
@@ -29,6 +29,7 @@
 #include "mlir/Analysis/Presburger/Matrix.h"
 #include "mlir/Analysis/Presburger/PresburgerRelation.h"
 #include "mlir/Analysis/Presburger/QuasiPolynomial.h"
+#include <bitset>
 #include <optional>
 
 namespace mlir {
@@ -91,6 +92,19 @@ GeneratingFunction unimodularConeGeneratingFunction(ParamPoint vertex, int sign,
 /// If there is no solution, return null.
 std::optional<ParamPoint> solveParametricEquations(FracMatrix equations);
 
+/// Given a list of possibly intersecting regions (PresburgerRelations) and
+/// the vertices active in each region, produce a pairwise disjoint list of
+/// regions (chambers) and identify the vertices active in each of these new
+/// regions.
+/// In the return type, the vertices are stored by their index in the
+/// `vertices` argument, i.e., the set {2, 3, 4} represents the vertex set
+/// {vertices[2], vertices[3], vertices[4]}.
+/// Note that here, by disjoint, we mean that the intersection is not
+/// full-dimensional.
+std::vector<std::pair<PresburgerRelation, std::vector<unsigned>>>
+chamberDecomposition(std::vector<PresburgerRelation> activeRegions,
+                     std::vector<ParamPoint> vertices);
+
 /// Compute the generating function corresponding to a polytope.
 /// All tangent cones of the polytope must be unimodular.
 std::vector<std::pair<PresburgerRelation, GeneratingFunction>>
diff --git a/mlir/include/mlir/Analysis/Presburger/Matrix.h b/mlir/include/mlir/Analysis/Presburger/Matrix.h
index 4a5bba17e204531..ae781a93a036ff6 100644
--- a/mlir/include/mlir/Analysis/Presburger/Matrix.h
+++ b/mlir/include/mlir/Analysis/Presburger/Matrix.h
@@ -20,6 +20,7 @@
 #include "llvm/ADT/ArrayRef.h"
 #include "llvm/Support/raw_ostream.h"
 
+#include <bitset>
 #include <cassert>
 
 namespace mlir {
@@ -190,6 +191,12 @@ class Matrix {
   Matrix<T> getSubMatrix(unsigned fromRow, unsigned toRow, unsigned fromColumn,
                          unsigned toColumn) const;
 
+  /// Split the rows of a matrix into two matrices according to which bits are
+  /// 1 and which are 0 in a given bitset.
+  /// The first matrix returned has the rows corresponding to 1 and the second
+  /// corresponding to 2.
+  std::pair<Matrix<T>, Matrix<T>> splitByBitset(std::bitset<16> indicator);
+
   /// Print the matrix.
   void print(raw_ostream &os) const;
   void dump() const;
@@ -303,6 +310,10 @@ class FracMatrix : public Matrix<Fraction> {
   // paper](https://www.cs.cmu.edu/~avrim/451f11/lectures/lect1129_LLL.pdf)
   // calls `y`, usually 3/4.
   void LLL(Fraction delta);
+
+  // Multiply each row of the matrix by the LCD of the denominators, thereby
+  // converting it to an integer matrix.
+  IntMatrix normalizeRows();
 };
 
 } // namespace presburger
diff --git a/mlir/lib/Analysis/Presburger/Barvinok.cpp b/mlir/lib/Analysis/Presburger/Barvinok.cpp
index 8c4de595b685c49..00a1ee97cd38de6 100644
--- a/mlir/lib/Analysis/Presburger/Barvinok.cpp
+++ b/mlir/lib/Analysis/Presburger/Barvinok.cpp
@@ -223,6 +223,92 @@ mlir::presburger::detail::solveParametricEquations(FracMatrix equations) {
   return vertex;
 }
 
+/// This is an implementation of the Clauss-Loechner algorithm for chamber
+/// decomposition.
+/// We maintain a list of pairwise disjoint chambers and their vertex-sets;
+/// we iterate over the vertex list, each time appending the vertex to the
+/// chambers where it is active and creating a new chamber if necessary.
+std::vector<std::pair<PresburgerRelation, std::vector<unsigned>>>
+mlir::presburger::detail::chamberDecomposition(
+    std::vector<PresburgerRelation> activeRegions,
+    std::vector<ParamPoint> vertices) {
+  // We maintain a list of regions and their associated vertex sets,
+  // initialized with the first vertex and its corresponding activity region.
+  std::vector<std::pair<PresburgerRelation, std::vector<unsigned>>> chambers = {
+      std::make_pair(activeRegions[0], std::vector({0u}))};
+  // Note that instead of storing lists of actual vertices, we store lists
+  // of indices. Thus the set {2, 3, 4} represents the vertex set
+  // {vertices[2], vertices[3], vertices[4]}.
+
+  std::vector<std::pair<PresburgerRelation, std::vector<unsigned>>> newChambers;
+
+  // We iterate over the vertex set.
+  // For each vertex v_j and its activity region R_j,
+  // we examine all the current chambers R_i.
+  // If R_j has a full-dimensional intersection with an existing chamber R_i,
+  // then that chamber is replaced by two new ones:
+  // 1. the intersection R_i \cap R_j, where v_j is active;
+  // 2. the difference R_i - R_j, where v_j is inactive.
+  // Once we have examined all R_i, we add a final chamber
+  // R_j - (union of all existing chambers),
+  // in which only v_j is active.
+  for (unsigned j = 1, e = vertices.size(); j < e; j++) {
+    newChambers.clear();
+
+    PresburgerRelation newRegion = activeRegions[j];
+    ParamPoint newVertex = vertices[j];
+
+    for (unsigned i = 0, f = chambers.size(); i < f; i++) {
+      auto [currentRegion, currentVertices] = chambers[i];
+
+      // First, we check if the intersection of R_j and R_i.
+      // It is a disjoint union of convex regions in the parameter space,
+      // and so we know that it is full-dimensional if any of the disjuncts
+      // is full-dimensional.
+      PresburgerRelation intersection = currentRegion.intersect(newRegion);
+      bool isFullDim = intersection.getNumRangeVars() == 0 ||
+                       llvm::any_of(intersection.getAllDisjuncts(),
+                                    [&](IntegerRelation disjunct) -> bool {
+                                      return disjunct.isFullDim();
+                                    });
+
+      // If the intersection is not full-dimensional, we do not modify
+      // the chamber list.
+      if (!isFullDim)
+        newChambers.push_back(chambers[i]);
+      else {
+        // If it is, we add the intersection and the difference as new chambers.
+        PresburgerRelation subtraction = currentRegion.subtract(newRegion);
+        newChambers.push_back(std::make_pair(subtraction, currentVertices));
+
+        currentVertices.push_back(j);
+        newChambers.push_back(std::make_pair(intersection, currentVertices));
+      }
+    }
+
+    // Finally we compute the chamber where only v_j is active by subtracting
+    // all existing chambers from R_j.
+    for (const std::pair<PresburgerRelation, std::vector<unsigned>> &chamber :
+         newChambers)
+      newRegion = newRegion.subtract(chamber.first);
+    newChambers.push_back(std::make_pair(newRegion, std::vector({j})));
+
+    // We filter `chambers` to remove empty regions.
+    chambers.clear();
+    for (const std::pair<PresburgerRelation, std::vector<unsigned>> &chamber :
+         newChambers) {
+      auto [r, v] = chamber;
+      bool isEmpty = llvm::all_of(
+          r.getAllDisjuncts(),
+          [&](IntegerRelation disjunct) -> bool { return disjunct.isEmpty(); });
+      if (!isEmpty)
+        chambers.push_back(chamber);
+    }
+  }
+
+  return chambers;
+}
+
 /// For a polytope expressed as a set of inequalities, compute the generating
 /// function corresponding to the number of lattice points present. This
 /// algorithm has three main steps:
@@ -282,25 +368,16 @@ mlir::presburger::detail::polytopeGeneratingFunction(PolyhedronH poly) {
     if (indicator.count() != numVars)
       continue;
 
-    // Collect the inequalities corresponding to the bits which are set.
-    IntMatrix subset(numVars, numVars + numSymbols + 1);
-    unsigned j1 = 0, j2 = 0;
-    for (unsigned i = 0; i < numIneqs; i++) {
-      if (indicator.test(i))
-        subset.setRow(j1++, poly.getInequality(i));
-
-      else {
-        // All other inequalities are stored in a2 and b2c2.
-        // These are column-wise splits of the inequalities;
-        // a2 stores the coefficients of the variables, and
-        // b2c2 stores the coefficients of the parameters and the constant term.
-        for (unsigned k = 0; k < numVars; k++)
-          a2(j2, k) = poly.atIneq(i, k);
-        for (unsigned k = numVars; k < numVars + numSymbols + 1; k++)
-          b2c2(j2, k - numVars) = poly.atIneq(i, k);
-        j2++;
-      }
-    }
+    // Collect the inequalities corresponding to the bits which are set
+    // and the remaining ones.
+    auto [subset, remainder] = poly.getInequalities().splitByBitset(indicator);
+    // All other inequalities are stored in a2 and b2c2.
+    // These are column-wise splits of the inequalities;
+    // a2 stores the coefficients of the variables, and
+    // b2c2 stores the coefficients of the parameters and the constant term.
+    a2 = FracMatrix(remainder.getSubMatrix(0, numIneqs - 1, 0, numVars - 1));
+    b2c2 = FracMatrix(
+        remainder.getSubMatrix(0, numIneqs - 1, numVars, numVars + numSymbols));
 
     // Find the vertex, if any, corresponding to the current subset of
     // inequalities.
@@ -337,22 +414,11 @@ mlir::presburger::detail::polytopeGeneratingFunction(PolyhedronH poly) {
 
     // We convert the representation of the active region to an integers-only
     // form so as to store it as an PresburgerRelation.
-    // We do this by taking the LCM of the denominators of all the coefficients
-    // and multiplying by it throughout.
-    IntMatrix activeRegionNorm = IntMatrix(numIneqs - numVars, numSymbols + 1);
+    IntMatrix activeRegionNorm = activeRegion.normalizeRows();
     IntegerRelation activeRegionRel =
         IntegerRelation(PresburgerSpace::getRelationSpace(0, numSymbols, 0, 0));
-    MPInt lcmDenoms = MPInt(1);
-    for (unsigned i = 0; i < numIneqs - numVars; i++) {
-      for (unsigned j = 0; j < numSymbols + 1; j++)
-        lcmDenoms = lcm(lcmDenoms, activeRegion(i, j).den);
-      for (unsigned j = 0; j < numSymbols + 1; j++)
-        activeRegionNorm(i, j) =
-            (activeRegion(i, j) * lcmDenoms).getAsInteger();
-
+    for (unsigned i = 0, e = activeRegion.getNumRows(); i < e; ++i)
       activeRegionRel.addInequality(activeRegionNorm.getRow(i));
-    }
-
     activeRegions.push_back(PresburgerRelation(activeRegionRel));
   }
 
@@ -361,80 +427,8 @@ mlir::presburger::detail::polytopeGeneratingFunction(PolyhedronH poly) {
   // property that no two of them have a full-dimensional intersection, i.e.,
   // they may share "faces" or "edges", but their intersection can only have
   // up to numVars-1 dimensions.
-
-  // We maintain a list of regions and their associated vertex sets,
-  // initialized with the first vertex and its corresponding activity region.
-  std::vector<std::pair<PresburgerRelation, std::vector<unsigned>>> chambers = {
-      std::make_pair(activeRegions[0], std::vector({0u}))};
-  // Note that instead of storing lists of actual vertices, we store lists
-  // of indices. Thus the set {2, 3, 4} represents the vertex set
-  // {vertices[2], vertices[3], vertices[4]}.
-
-  std::vector<std::pair<PresburgerRelation, std::vector<unsigned>>> newChambers;
-
-  // We iterate over the vertex set.
-  // For each vertex v_j and its activity region R_j,
-  // we examine all the current chambers R_i.
-  // If R_j has a full-dimensional intersection with an existing chamber R_i,
-  // then that chamber is replaced by two new ones:
-  // 1. the intersection R_i \cap R_j, where v_j is active;
-  // 2. the difference R_i - R_j, where v_j is inactive.
-  // Once we have examined all R_i, we add a final chamber
-  // R_j - (union of all existing chambers),
-  // in which only v_j is active.
-  for (unsigned j = 1, e = vertices.size(); j < e; j++) {
-    newChambers.clear();
-
-    PresburgerRelation newRegion = activeRegions[j];
-    ParamPoint newVertex = vertices[j];
-
-    for (unsigned i = 0, f = chambers.size(); i < f; i++) {
-      auto [currentRegion, currentVertices] = chambers[i];
-
-      // First, we check if the intersection of R_j and R_i.
-      // It is a disjoint union of convex regions in the parameter space,
-      // and so we know that it is full-dimensional if any of the disjuncts
-      // is full-dimensional.
-      PresburgerRelation intersection = currentRegion.intersect(newRegion);
-      bool isFullDim = numSymbols == 0 ||
-                       llvm::any_of(intersection.getAllDisjuncts(),
-                                    [&](IntegerRelation disjunct) -> bool {
-                                      return disjunct.isFullDim();
-                                    });
-
-      // If the intersection is not full-dimensional, we do not modify
-      // the chamber list.
-      if (!isFullDim)
-        newChambers.push_back(chambers[i]);
-      else {
-        // If it is, we add the intersection and the difference as new chambers.
-        PresburgerRelation subtraction = currentRegion.subtract(newRegion);
-        newChambers.push_back(std::make_pair(subtraction, currentVertices));
-
-        currentVertices.push_back(j);
-        newChambers.push_back(std::make_pair(intersection, currentVertices));
-      }
-    }
-
-    // Finally we compute the chamber where only v_j is active by subtracting
-    // all existing chambers from R_j.
-    for (const std::pair<PresburgerRelation, std::vector<unsigned>> &chamber :
-         newChambers)
-      newRegion = newRegion.subtract(chamber.first);
-    newChambers.push_back(std::make_pair(newRegion, std::vector({j})));
-
-    // We filter `chambers` to remove empty regions.
-    chambers.clear();
-    for (const std::pair<PresburgerRelation, std::vector<unsigned>> &chamber :
-         newChambers) {
-      auto [r, v] = chamber;
-      bool isEmpty = llvm::all_of(
-          r.getAllDisjuncts(),
-          [&](IntegerRelation disjunct) -> bool { return disjunct.isEmpty(); });
-      if (!isEmpty)
-        chambers.push_back(chamber);
-    }
-  }
+  std::vector<std::pair<PresburgerRelation, std::vector<unsigned>>> chambers =
+      chamberDecomposition(activeRegions, vertices);
 
   // Now, we compute the generating function. For each chamber, we iterate over
   // the vertices active in it, and compute the generating function for each
diff --git a/mlir/lib/Analysis/Presburger/Matrix.cpp b/mlir/lib/Analysis/Presburger/Matrix.cpp
index be8dac636e0af61..099cb6fa09e767f 100644
--- a/mlir/lib/Analysis/Presburger/Matrix.cpp
+++ b/mlir/lib/Analysis/Presburger/Matrix.cpp
@@ -376,6 +376,21 @@ void Matrix<T>::print(raw_ostream &os) const {
   }
 }
 
+/// We iterate over the `indicator` bitset, checking each bit. If a bit is 1,
+/// we append it to one matrix, and if it is zero, we append it to the other.
+template <typename T>
+std::pair<Matrix<T>, Matrix<T>>
+Matrix<T>::splitByBitset(std::bitset<16> indicator) {
+  Matrix<T> rowsForOne(0, nColumns), rowsForZero(0, nColumns);
+  for (unsigned i = 0; i < nRows; i++) {
+    if (indicator.test(i))
+      rowsForOne.appendExtraRow(getRow(i));
+    else
+      rowsForZero.appendExtraRow(getRow(i));
+  }
+  return std::make_pair(rowsForOne, rowsForZero);
+}
+
 template <typename T>
 void Matrix<T>::dump() const {
   print(llvm::errs());
@@ -710,3 +725,20 @@ void FracMatrix::LLL(Fraction delta) {
     }
   }
 }
+
+IntMatrix FracMatrix::normalizeRows() {
+  unsigned numRows = getNumRows();
+  unsigned numColumns = getNumColumns();
+  IntMatrix normalized(numRows, numColumns);
+
+  MPInt lcmDenoms = MPInt(1);
+  for (unsigned i = 0; i < numRows; i++) {
+    // For a row, first compute the LCM of the denominators.
+    for (unsigned j = 0; j < numColumns; j++)
+      lcmDenoms = lcm(lcmDenoms, at(i, j).den);
+    // Then, multiply by it throughout and convert to integers.
+    for (unsigned j = 0; j < numColumns; j++)
+      normalized(i, j) = (at(i, j) * lcmDenoms).getAsInteger();
+  }
+  return normalized;
+}

>From 6f33bc045e578155dc538877d7fa19a8e859c475 Mon Sep 17 00:00:00 2001
From: Abhinav271828 <abhinav.m at research.iiit.ac.in>
Date: Tue, 23 Jan 2024 18:21:39 +0530
Subject: [PATCH 14/26] Fixes

---
 mlir/include/mlir/Analysis/Presburger/Barvinok.h |  4 ++--
 mlir/lib/Analysis/Presburger/Barvinok.cpp        | 10 +++++-----
 mlir/lib/Analysis/Presburger/Matrix.cpp          |  2 +-
 3 files changed, 8 insertions(+), 8 deletions(-)

diff --git a/mlir/include/mlir/Analysis/Presburger/Barvinok.h b/mlir/include/mlir/Analysis/Presburger/Barvinok.h
index 34055e93242c4ca..558498e99104095 100644
--- a/mlir/include/mlir/Analysis/Presburger/Barvinok.h
+++ b/mlir/include/mlir/Analysis/Presburger/Barvinok.h
@@ -49,7 +49,7 @@ using PolyhedronV = IntMatrix;
 using ConeH = PolyhedronH;
 using ConeV = PolyhedronV;
 
-inline ConeH defineHRep(int numVars) {
+inline ConeH defineHRep(int numVars, int numSymbols = 0) {
   // We don't distinguish between domain and range variables, so
   // we set the number of domain variables as 0 and the number of
   // range variables as the number of actual variables.
@@ -57,7 +57,7 @@ inline ConeH defineHRep(int numVars) {
   // (existentially quantified) variables.
   // Once the cone is defined, we use `addInequality()` to set inequalities.
   return ConeH(PresburgerSpace::getSetSpace(/*numDims=*/numVars,
-                                            /*numSymbols=*/0,
+                                            /*numSymbols=*/numSymbols,
                                             /*numLocals=*/0));
 }
 
diff --git a/mlir/lib/Analysis/Presburger/Barvinok.cpp b/mlir/lib/Analysis/Presburger/Barvinok.cpp
index 00a1ee97cd38de6..4b32489943bb0cf 100644
--- a/mlir/lib/Analysis/Presburger/Barvinok.cpp
+++ b/mlir/lib/Analysis/Presburger/Barvinok.cpp
@@ -360,8 +360,7 @@ mlir::presburger::detail::polytopeGeneratingFunction(PolyhedronH poly) {
   // those numbers whose binary representation has numVars 1's and the
   // rest 0's.
   unsigned upperBound = ((1ul << numVars) - 1ul) << (numIneqs - numVars);
-  for (std::bitset<16> indicator(((1ul << numVars) - 1ul)
-                                 << (numIneqs - numVars));
+  for (std::bitset<16> indicator(upperBound);
        indicator.to_ulong() <= upperBound;
        indicator = std::bitset<16>(indicator.to_ulong() - 1)) {
 
@@ -375,9 +374,10 @@ mlir::presburger::detail::polytopeGeneratingFunction(PolyhedronH poly) {
     // These are column-wise splits of the inequalities;
     // a2 stores the coefficients of the variables, and
     // b2c2 stores the coefficients of the parameters and the constant term.
-    a2 = FracMatrix(remainder.getSubMatrix(0, numIneqs - 1, 0, numVars - 1));
-    b2c2 = FracMatrix(
-        remainder.getSubMatrix(0, numIneqs - 1, numVars, numVars + numSymbols));
+    a2 = FracMatrix(
+        remainder.getSubMatrix(0, numIneqs - numVars - 1, 0, numVars - 1));
+    b2c2 = FracMatrix(remainder.getSubMatrix(0, numIneqs - numVars - 1, numVars,
+                                             numVars + numSymbols));
 
     // Find the vertex, if any, corresponding to the current subset of
     // inequalities.
diff --git a/mlir/lib/Analysis/Presburger/Matrix.cpp b/mlir/lib/Analysis/Presburger/Matrix.cpp
index 099cb6fa09e767f..35af156552c34d9 100644
--- a/mlir/lib/Analysis/Presburger/Matrix.cpp
+++ b/mlir/lib/Analysis/Presburger/Matrix.cpp
@@ -363,7 +363,7 @@ Matrix<T> Matrix<T>::getSubMatrix(unsigned fromRow, unsigned toRow,
   Matrix<T> subMatrix(toRow - fromRow + 1, toColumn - fromColumn + 1);
   for (unsigned i = fromRow; i <= toRow; i++)
     for (unsigned j = fromColumn; j <= toColumn; j++)
-      subMatrix(i, j) = at(i, j);
+      subMatrix(i - fromRow, j - fromColumn) = at(i, j);
   return subMatrix;
 }
 

>From 93ed01375b427732ea64f822e7b737d7bb81369e Mon Sep 17 00:00:00 2001
From: Abhinav271828 <abhinav.m at research.iiit.ac.in>
Date: Tue, 23 Jan 2024 18:45:02 +0530
Subject: [PATCH 15/26] Fix index

---
 mlir/lib/Analysis/Presburger/Barvinok.cpp     |  2 +-
 .../Analysis/Presburger/BarvinokTest.cpp      | 20 +++++++++++++++++--
 2 files changed, 19 insertions(+), 3 deletions(-)

diff --git a/mlir/lib/Analysis/Presburger/Barvinok.cpp b/mlir/lib/Analysis/Presburger/Barvinok.cpp
index 4b32489943bb0cf..82fbee5a05b44e9 100644
--- a/mlir/lib/Analysis/Presburger/Barvinok.cpp
+++ b/mlir/lib/Analysis/Presburger/Barvinok.cpp
@@ -216,7 +216,7 @@ mlir::presburger::detail::solveParametricEquations(FracMatrix equations) {
   // We copy these columns and return them.
   ParamPoint vertex =
       equations.getSubMatrix(/*fromRow=*/0, /*toRow=*/numEqs - 1,
-                             /*fromRow=*/numEqs, /*toRow=*/numCols);
+                             /*fromColumn=*/numEqs, /*toColumn=*/numCols - 1);
   for (unsigned i = 0; i < numEqs; ++i)
     vertex.negateRow(i);
 
diff --git a/mlir/unittests/Analysis/Presburger/BarvinokTest.cpp b/mlir/unittests/Analysis/Presburger/BarvinokTest.cpp
index 919aaa7a4285930..9770d03c90a9f0e 100644
--- a/mlir/unittests/Analysis/Presburger/BarvinokTest.cpp
+++ b/mlir/unittests/Analysis/Presburger/BarvinokTest.cpp
@@ -125,7 +125,7 @@ TEST(BarvinokTest, getCoefficientInRationalFunction) {
   EXPECT_EQ(coeff.getConstantTerm(), Fraction(55, 64));
 }
 
-TEST(BarvinokTest, computeNumTerms) {
+TEST(BarvinokTest, computeNumTermsCone) {
   // The following test is taken from
   // Verdoolaege, Sven, et al. "Counting integer points in parametric
   // polytopes using Barvinok's rational functions." Algorithmica 48 (2007):
@@ -233,4 +233,20 @@ TEST(BarvinokTest, computeNumTerms) {
     for (unsigned j = 0; j < 2; j++)
       for (unsigned k = 0; k < 2; k++)
         EXPECT_EQ(count[i][j][k], 1);
-}
\ No newline at end of file
+}
+
+TEST(BarvinokTest, computeNumTermsPolytope) {
+  IntMatrix ineqs = makeIntMatrix(6, 4, {{1, 0, 0, 0},
+                                         {0, 1, 0, 0},
+                                         {0, 0, 1, 0},
+                                         {-1, 0, 0, 1},
+                                         {0, -1, 0, 1},
+                                         {0, 0, -1, 1}});
+    PolyhedronH poly = defineHRep(3);
+    for (unsigned i = 0; i < 6; i++)
+        poly.addInequality(ineqs.getRow(i));
+
+  std::vector<std::pair<PresburgerRelation, GeneratingFunction>> count = polytopeGeneratingFunction(poly);
+  EXPECT_EQ(count.size(), 1u);
+
+}

>From fe73f449fd801c140902481515bd1e6a0c825fca Mon Sep 17 00:00:00 2001
From: Abhinav271828 <abhinav.m at research.iiit.ac.in>
Date: Tue, 23 Jan 2024 18:48:49 +0530
Subject: [PATCH 16/26] Minor fix

---
 mlir/lib/Analysis/Presburger/Barvinok.cpp | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/mlir/lib/Analysis/Presburger/Barvinok.cpp b/mlir/lib/Analysis/Presburger/Barvinok.cpp
index 82fbee5a05b44e9..ae66bedd21cd2f6 100644
--- a/mlir/lib/Analysis/Presburger/Barvinok.cpp
+++ b/mlir/lib/Analysis/Presburger/Barvinok.cpp
@@ -85,7 +85,7 @@ GeneratingFunction mlir::presburger::detail::unimodularConeGeneratingFunction(
   //                                       [-1 -1/2 1]
 
   // `cone` must be unimodular.
-  assert(getIndex(getDual(cone)) == 1 && "input cone is not unimodular!");
+  assert(abs(getIndex(getDual(cone))) == 1 && "input cone is not unimodular!");
 
   unsigned numVar = cone.getNumVars();
   unsigned numIneq = cone.getNumInequalities();
@@ -443,11 +443,12 @@ mlir::presburger::detail::polytopeGeneratingFunction(PolyhedronH poly) {
       // We collect the inequalities corresponding to each vertex.
       // We only need the coefficients of the variables (NOT the parameters)
       // as the generating function only depends on these.
+      // We translate the cones to be pointed at the origin by making the
+      // constant terms zero.
       ConeH tgtCone = defineHRep(numVars);
       for (unsigned j = 0; j < numVars; j++) {
         for (unsigned k = 0; k < numVars; k++)
           ineq[k] = subsets[i](j, k);
-        ineq[numVars] = subsets[i](j, numVars + numSymbols);
         tgtCone.addInequality(ineq);
       }
       // We assume that the tangent cone is unimodular.

>From 2d8512f8cf5bb576f33380289051a97afd544b5c Mon Sep 17 00:00:00 2001
From: Abhinav271828 <abhinav.m at research.iiit.ac.in>
Date: Tue, 23 Jan 2024 19:32:06 +0530
Subject: [PATCH 17/26] Add test for pgf

---
 .../Analysis/Presburger/BarvinokTest.cpp      | 48 ++++++++++++++-----
 1 file changed, 37 insertions(+), 11 deletions(-)

diff --git a/mlir/unittests/Analysis/Presburger/BarvinokTest.cpp b/mlir/unittests/Analysis/Presburger/BarvinokTest.cpp
index 9770d03c90a9f0e..cd616fb7bbf3cf5 100644
--- a/mlir/unittests/Analysis/Presburger/BarvinokTest.cpp
+++ b/mlir/unittests/Analysis/Presburger/BarvinokTest.cpp
@@ -235,18 +235,44 @@ TEST(BarvinokTest, computeNumTermsCone) {
         EXPECT_EQ(count[i][j][k], 1);
 }
 
+/// We define some simple polyhedra with unimodular tangent cones and verify
+/// that the returned generating functions correspond to those calculated by
+/// hand.
 TEST(BarvinokTest, computeNumTermsPolytope) {
-  IntMatrix ineqs = makeIntMatrix(6, 4, {{1, 0, 0, 0},
-                                         {0, 1, 0, 0},
-                                         {0, 0, 1, 0},
-                                         {-1, 0, 0, 1},
-                                         {0, -1, 0, 1},
-                                         {0, 0, -1, 1}});
-    PolyhedronH poly = defineHRep(3);
-    for (unsigned i = 0; i < 6; i++)
-        poly.addInequality(ineqs.getRow(i));
-
-  std::vector<std::pair<PresburgerRelation, GeneratingFunction>> count = polytopeGeneratingFunction(poly);
+  IntMatrix ineqs = makeIntMatrix(6, 4,
+                                  {{1, 0, 0, 0},
+                                   {0, 1, 0, 0},
+                                   {0, 0, 1, 0},
+                                   {-1, 0, 0, 1},
+                                   {0, -1, 0, 1},
+                                   {0, 0, -1, 1}});
+  PolyhedronH poly = defineHRep(3);
+  for (unsigned i = 0; i < 6; i++)
+    poly.addInequality(ineqs.getRow(i));
+
+  std::vector<std::pair<PresburgerRelation, GeneratingFunction>> count =
+      polytopeGeneratingFunction(poly);
+  // There is only one chamber, as it is non-parametric.
   EXPECT_EQ(count.size(), 1u);
 
+  GeneratingFunction gf = count[0].second;
+  EXPECT_EQ_REPR_GENERATINGFUNCTION(
+      gf,
+      GeneratingFunction(
+          0, {1, 1, 1, 1, 1, 1, 1, 1},
+          {makeFracMatrix(1, 3, {{1, 1, 1}}), makeFracMatrix(1, 3, {{0, 1, 1}}),
+           makeFracMatrix(1, 3, {{0, 1, 1}}), makeFracMatrix(1, 3, {{0, 0, 1}}),
+           makeFracMatrix(1, 3, {{0, 1, 1}}), makeFracMatrix(1, 3, {{0, 0, 1}}),
+           makeFracMatrix(1, 3, {{0, 0, 1}}),
+           makeFracMatrix(1, 3, {{0, 0, 0}})},
+          {{{-1, 0, 0}, {0, -1, 0}, {0, 0, -1}},
+           {{1, 0, 0}, {0, -1, 0}, {0, 0, -1}},
+           {{0, 1, 0}, {-1, 0, 0}, {0, 0, -1}},
+           {{1, 0, 0}, {0, 1, 0}, {0, 0, -1}},
+           {{0, 0, 1}, {-1, 0, 0}, {0, -1, 0}},
+           {{1, 0, 0}, {0, 0, 1}, {0, -1, 0}},
+           {{0, 1, 0}, {0, 0, 1}, {-1, 0, 0}},
+           {{1, 0, 0}, {0, 1, 0}, {0, 0, 1}}}));
+
+
 }

>From 7c2906dcb6150f78a89188d54bf5c9c7d7d7a9d2 Mon Sep 17 00:00:00 2001
From: Abhinav271828 <abhinav.m at research.iiit.ac.in>
Date: Tue, 23 Jan 2024 19:34:58 +0530
Subject: [PATCH 18/26] Fix doc

---
 mlir/include/mlir/Analysis/Presburger/Barvinok.h | 10 ++++++----
 1 file changed, 6 insertions(+), 4 deletions(-)

diff --git a/mlir/include/mlir/Analysis/Presburger/Barvinok.h b/mlir/include/mlir/Analysis/Presburger/Barvinok.h
index 558498e99104095..e6d4af510b1ca68 100644
--- a/mlir/include/mlir/Analysis/Presburger/Barvinok.h
+++ b/mlir/include/mlir/Analysis/Presburger/Barvinok.h
@@ -49,14 +49,16 @@ using PolyhedronV = IntMatrix;
 using ConeH = PolyhedronH;
 using ConeV = PolyhedronV;
 
-inline ConeH defineHRep(int numVars, int numSymbols = 0) {
+inline PolyhedronH defineHRep(int numVars, int numSymbols = 0) {
   // We don't distinguish between domain and range variables, so
   // we set the number of domain variables as 0 and the number of
   // range variables as the number of actual variables.
-  // There are no symbols (we don't work with parametric cones) and no local
-  // (existentially quantified) variables.
+  // numSymbols is the number of parameters; it is
+  // There are no local (existentially quantified) variables.
+  // The number of symbols is the number of parameters. By default, we consider
+  // nonparametric polyhedra.
   // Once the cone is defined, we use `addInequality()` to set inequalities.
-  return ConeH(PresburgerSpace::getSetSpace(/*numDims=*/numVars,
+  return PolyhedronH(PresburgerSpace::getSetSpace(/*numDims=*/numVars,
                                             /*numSymbols=*/numSymbols,
                                             /*numLocals=*/0));
 }

>From 8b1169659263ca186f8dfc2b68ad5ef9f343caad Mon Sep 17 00:00:00 2001
From: Abhinav271828 <abhinav.m at research.iiit.ac.in>
Date: Tue, 23 Jan 2024 19:58:04 +0530
Subject: [PATCH 19/26] Fix removeTrivialEquality()

---
 .../Analysis/Presburger/IntegerRelation.cpp   | 29 ++++++++++---------
 1 file changed, 15 insertions(+), 14 deletions(-)

diff --git a/mlir/lib/Analysis/Presburger/IntegerRelation.cpp b/mlir/lib/Analysis/Presburger/IntegerRelation.cpp
index 1683f5cf41aee2e..5895b927435a424 100644
--- a/mlir/lib/Analysis/Presburger/IntegerRelation.cpp
+++ b/mlir/lib/Analysis/Presburger/IntegerRelation.cpp
@@ -2499,30 +2499,31 @@ void IntegerRelation::printSpace(raw_ostream &os) const {
 }
 
 void IntegerRelation::removeTrivialEqualities() {
-  bool flag;
-  for (unsigned i = 0, e = getNumInequalities(); i < e; i++) {
-    flag = true;
-    for (unsigned j = 0, f = getNumVars(); j < f + 1; j++)
-      if (atEq(i, j) != 0)
-        flag = false;
-    if (flag)
-      removeEquality(i);
+  std::vector<bool> isTrivialEquality;
+  for (int i = 0, e = getNumEqualities(); i < e; ++i) {
+    bool currentIsTrivial =
+        llvm::all_of(getEquality(i), [&](MPInt n) -> bool { return (n == 0); });
+    isTrivialEquality.push_back(currentIsTrivial);
   }
-}
 
-bool IntegerRelation::isFullDim() {
-  removeTrivialEqualities();
+  unsigned pos = 0;
+  for (unsigned r = 0, e = getNumEqualities(); r < e; r++)
+    if (!isTrivialEquality[r])
+      equalities.copyRow(r, pos++);
 
-  unsigned e = getNumInequalities();
+  equalities.resizeVertically(pos);
+}
 
+bool IntegerRelation::isFullDim() {
   // If there is a non-trivial equality, the space cannot be full-dimensional.
-  if (e > 0)
+  removeTrivialEqualities();
+  if (getNumEqualities() > 0)
     return false;
 
   // If along the direction of any of the inequalities, the upper and lower
   // optima are the same, then the region is not full-dimensional.
   Simplex simplex(*this);
-  for (unsigned i = 0; i < e; i++) {
+  for (unsigned i = 0, e = getNumInequalities(); i < e; i++) {
     auto ineq = inequalities.getRow(i);
     auto upOpt = simplex.computeOptimum(Simplex::Direction::Up, ineq);
     auto downOpt = simplex.computeOptimum(Simplex::Direction::Down, ineq);

>From c9448e38f24f92a3752d01ce7715f1ff87229fe0 Mon Sep 17 00:00:00 2001
From: Abhinav271828 <abhinav.m at research.iiit.ac.in>
Date: Tue, 23 Jan 2024 20:03:52 +0530
Subject: [PATCH 20/26] Formatting

---
 mlir/include/mlir/Analysis/Presburger/Barvinok.h    |  4 ++--
 mlir/unittests/Analysis/Presburger/BarvinokTest.cpp | 10 ++++++++++
 2 files changed, 12 insertions(+), 2 deletions(-)

diff --git a/mlir/include/mlir/Analysis/Presburger/Barvinok.h b/mlir/include/mlir/Analysis/Presburger/Barvinok.h
index e6d4af510b1ca68..3a80d201730f5df 100644
--- a/mlir/include/mlir/Analysis/Presburger/Barvinok.h
+++ b/mlir/include/mlir/Analysis/Presburger/Barvinok.h
@@ -59,8 +59,8 @@ inline PolyhedronH defineHRep(int numVars, int numSymbols = 0) {
   // nonparametric polyhedra.
   // Once the cone is defined, we use `addInequality()` to set inequalities.
   return PolyhedronH(PresburgerSpace::getSetSpace(/*numDims=*/numVars,
-                                            /*numSymbols=*/numSymbols,
-                                            /*numLocals=*/0));
+                                                  /*numSymbols=*/numSymbols,
+                                                  /*numLocals=*/0));
 }
 
 /// Get the index of a cone, i.e., the volume of the parallelepiped
diff --git a/mlir/unittests/Analysis/Presburger/BarvinokTest.cpp b/mlir/unittests/Analysis/Presburger/BarvinokTest.cpp
index cd616fb7bbf3cf5..81c9b79bb14a338 100644
--- a/mlir/unittests/Analysis/Presburger/BarvinokTest.cpp
+++ b/mlir/unittests/Analysis/Presburger/BarvinokTest.cpp
@@ -239,6 +239,7 @@ TEST(BarvinokTest, computeNumTermsCone) {
 /// that the returned generating functions correspond to those calculated by
 /// hand.
 TEST(BarvinokTest, computeNumTermsPolytope) {
+  // A cube of side 1.
   IntMatrix ineqs = makeIntMatrix(6, 4,
                                   {{1, 0, 0, 0},
                                    {0, 1, 0, 0},
@@ -274,5 +275,14 @@ TEST(BarvinokTest, computeNumTermsPolytope) {
            {{0, 1, 0}, {0, 0, 1}, {-1, 0, 0}},
            {{1, 0, 0}, {0, 1, 0}, {0, 0, 1}}}));
 
+  // A right-angled triangle with side p.
+  ineqs = makeIntMatrix(3, 4, {{1, 0, 0, 0}, {0, 1, 0, 0}, {-1, -1, 1, 0}});
+  poly = defineHRep(2, 1);
+  for (unsigned i = 0; i < 3; i++)
+    poly.addInequality(ineqs.getRow(i));
+  
+  count = polytopeGeneratingFunction(poly);
+  // There is only one chamber: p ≥ 0
+  EXPECT_EQ(count.size(), 1u);
 
 }

>From 8a21ff3937ff0597514047318d0585e1225148de Mon Sep 17 00:00:00 2001
From: Abhinav271828 <abhinav.m at research.iiit.ac.in>
Date: Tue, 23 Jan 2024 20:25:30 +0530
Subject: [PATCH 21/26] Add test for pgf

---
 mlir/unittests/Analysis/Presburger/BarvinokTest.cpp | 10 +++++++++-
 1 file changed, 9 insertions(+), 1 deletion(-)

diff --git a/mlir/unittests/Analysis/Presburger/BarvinokTest.cpp b/mlir/unittests/Analysis/Presburger/BarvinokTest.cpp
index 81c9b79bb14a338..09301d5e54e0e0c 100644
--- a/mlir/unittests/Analysis/Presburger/BarvinokTest.cpp
+++ b/mlir/unittests/Analysis/Presburger/BarvinokTest.cpp
@@ -280,9 +280,17 @@ TEST(BarvinokTest, computeNumTermsPolytope) {
   poly = defineHRep(2, 1);
   for (unsigned i = 0; i < 3; i++)
     poly.addInequality(ineqs.getRow(i));
-  
+
   count = polytopeGeneratingFunction(poly);
   // There is only one chamber: p ≥ 0
   EXPECT_EQ(count.size(), 1u);
 
+  gf = count[0].second;
+  EXPECT_EQ_REPR_GENERATINGFUNCTION(
+      gf, GeneratingFunction(
+              1, {1, 1, 1},
+              {makeFracMatrix(2, 2, {{0, 1}, {0, 0}}),
+               makeFracMatrix(2, 2, {{0, 1}, {0, 0}}),
+               makeFracMatrix(2, 2, {{0, 0}, {0, 0}})},
+              {{{-1, 1}, {-1, 0}}, {{1, -1}, {0, -1}}, {{1, 0}, {0, 1}}}));
 }

>From 57f4da54904e51850dff11e416abd13096b4baa6 Mon Sep 17 00:00:00 2001
From: Abhinav271828 <abhinav.m at research.iiit.ac.in>
Date: Wed, 24 Jan 2024 03:04:16 +0530
Subject: [PATCH 22/26] Minor fixes

---
 .../mlir/Analysis/Presburger/Barvinok.h       |  4 +-
 .../Analysis/Presburger/GeneratingFunction.h  |  2 +-
 .../Analysis/Presburger/IntegerRelation.h     |  2 +
 .../mlir/Analysis/Presburger/Simplex.h        |  7 +--
 mlir/lib/Analysis/Presburger/Barvinok.cpp     | 39 ++++++++--------
 .../Analysis/Presburger/IntegerRelation.cpp   | 45 ++++++-------------
 mlir/lib/Analysis/Presburger/Matrix.cpp       |  4 +-
 mlir/lib/Analysis/Presburger/Simplex.cpp      | 18 +++++---
 8 files changed, 58 insertions(+), 63 deletions(-)

diff --git a/mlir/include/mlir/Analysis/Presburger/Barvinok.h b/mlir/include/mlir/Analysis/Presburger/Barvinok.h
index 3a80d201730f5df..4d3f643187247f2 100644
--- a/mlir/include/mlir/Analysis/Presburger/Barvinok.h
+++ b/mlir/include/mlir/Analysis/Presburger/Barvinok.h
@@ -104,8 +104,8 @@ std::optional<ParamPoint> solveParametricEquations(FracMatrix equations);
 /// Note that here, by disjoint, we mean that the intersection is not
 /// full-dimensional.
 std::vector<std::pair<PresburgerRelation, std::vector<unsigned>>>
-chamberDecomposition(std::vector<PresburgerRelation> activeRegions,
-                     std::vector<ParamPoint> vertices);
+computeChamberDecomposition(std::vector<PresburgerRelation> activeRegions,
+                            std::vector<ParamPoint> vertices);
 
 /// Compute the generating function corresponding to a polytope.
 /// All tangent cones of the polytope must be unimodular.
diff --git a/mlir/include/mlir/Analysis/Presburger/GeneratingFunction.h b/mlir/include/mlir/Analysis/Presburger/GeneratingFunction.h
index 35f15319b4568eb..db5b6b6a9591866 100644
--- a/mlir/include/mlir/Analysis/Presburger/GeneratingFunction.h
+++ b/mlir/include/mlir/Analysis/Presburger/GeneratingFunction.h
@@ -72,7 +72,7 @@ class GeneratingFunction {
     return denominators;
   }
 
-  GeneratingFunction operator+(GeneratingFunction gf) const {
+  GeneratingFunction operator+(const GeneratingFunction &gf) const {
     assert(numParam == gf.getNumParams() &&
            "two generating functions with different numbers of parameters "
            "cannot be added!");
diff --git a/mlir/include/mlir/Analysis/Presburger/IntegerRelation.h b/mlir/include/mlir/Analysis/Presburger/IntegerRelation.h
index b60d68aab4cf405..3e76e4bd2cc0ce2 100644
--- a/mlir/include/mlir/Analysis/Presburger/IntegerRelation.h
+++ b/mlir/include/mlir/Analysis/Presburger/IntegerRelation.h
@@ -716,6 +716,8 @@ class IntegerRelation {
 
   // Verify whether the relation is full-dimensional, i.e.,
   // no equality holds for the relation.
+  // If there are no variables, it always returns true.
+  // If there is at least one variable and the space is empty, it returns false.
   bool isFullDim();
 
   void print(raw_ostream &os) const;
diff --git a/mlir/include/mlir/Analysis/Presburger/Simplex.h b/mlir/include/mlir/Analysis/Presburger/Simplex.h
index 1939a41e52bc6d1..8be9942a3dab922 100644
--- a/mlir/include/mlir/Analysis/Presburger/Simplex.h
+++ b/mlir/include/mlir/Analysis/Presburger/Simplex.h
@@ -771,9 +771,10 @@ class Simplex : public SimplexBase {
   std::pair<MaybeOptimum<MPInt>, MaybeOptimum<MPInt>>
   computeIntegerBounds(ArrayRef<MPInt> coeffs);
 
-  /// Returns false if the given equality has distinct upper and lower bounds
-  /// in the simplex.
-  bool isValidEquality(const ArrayRef<MPInt> coeffs);
+  /// Check if the simplex takes only one rational value along the
+  /// direction of `coeffs`.
+  /// `this` must be nonempty.
+  bool isFlatAlong(const ArrayRef<MPInt> coeffs);
 
   /// Returns true if the polytope is unbounded, i.e., extends to infinity in
   /// some direction. Otherwise, returns false.
diff --git a/mlir/lib/Analysis/Presburger/Barvinok.cpp b/mlir/lib/Analysis/Presburger/Barvinok.cpp
index ae66bedd21cd2f6..70ec9da4d530f3d 100644
--- a/mlir/lib/Analysis/Presburger/Barvinok.cpp
+++ b/mlir/lib/Analysis/Presburger/Barvinok.cpp
@@ -161,33 +161,33 @@ std::optional<ParamPoint>
 mlir::presburger::detail::solveParametricEquations(FracMatrix equations) {
   // equations is a d x (d + p + 1) matrix.
   // Each row represents an equation.
-  unsigned numEqs = equations.getNumRows();
+  unsigned d = equations.getNumRows();
   unsigned numCols = equations.getNumColumns();
 
   // If the determinant is zero, there is no unique solution.
   // Thus we return null.
-  if (FracMatrix(equations.getSubMatrix(/*fromRow=*/0, /*toRow=*/numEqs - 1,
+  if (FracMatrix(equations.getSubMatrix(/*fromRow=*/0, /*toRow=*/d - 1,
                                         /*fromColumn=*/0,
-                                        /*toColumn=*/numEqs - 1))
+                                        /*toColumn=*/d - 1))
           .determinant() == 0)
     return std::nullopt;
 
-  for (unsigned i = 0; i < numEqs; ++i) {
+  for (unsigned i = 0; i < d; ++i) {
     // First ensure that the diagonal element is nonzero, by swapping
     // it with a row that is non-zero at column i.
-    if (equations(i, i) == 0) {
-      for (unsigned j = i + 1; j < numEqs; ++j) {
-        if (equations(j, i) == 0)
-          continue;
-        equations.swapRows(j, i);
-        break;
-      }
+    if (equations(i, i) != 0)
+      continue;
+    for (unsigned j = i + 1; j < d; ++j) {
+      if (equations(j, i) == 0)
+        continue;
+      equations.swapRows(j, i);
+      break;
     }
 
     Fraction diagElement = equations(i, i);
 
     // Apply row operations to make all elements except the diagonal to zero.
-    for (unsigned j = 0; j < numEqs; ++j) {
+    for (unsigned j = 0; j < d; ++j) {
       if (i == j)
         continue;
       if (equations(j, i) == 0)
@@ -195,12 +195,13 @@ mlir::presburger::detail::solveParametricEquations(FracMatrix equations) {
       // Apply row operations to make element (j, i) zero by subtracting the
       // ith row, appropriately scaled.
       Fraction currentElement = equations(j, i);
-      equations.addToRow(j, equations.getRow(i), -currentElement / diagElement);
+      equations.addToRow(/*sourceRow=*/i, /*targetRow=*/j,
+                         -currentElement / diagElement);
     }
   }
 
   // Rescale diagonal elements to 1.
-  for (unsigned i = 0; i < numEqs; ++i) {
+  for (unsigned i = 0; i < d; ++i) {
     Fraction diagElement = equations(i, i);
     for (unsigned j = 0; j < numCols; ++j)
       equations(i, j) = equations(i, j) / diagElement;
@@ -215,9 +216,9 @@ mlir::presburger::detail::solveParametricEquations(FracMatrix equations) {
   // and so we return the negation of the last p + 1 columns of the matrix.
   // We copy these columns and return them.
   ParamPoint vertex =
-      equations.getSubMatrix(/*fromRow=*/0, /*toRow=*/numEqs - 1,
-                             /*fromColumn=*/numEqs, /*toColumn=*/numCols - 1);
-  for (unsigned i = 0; i < numEqs; ++i)
+      equations.getSubMatrix(/*fromRow=*/0, /*toRow=*/d - 1,
+                             /*fromColumn=*/d, /*toColumn=*/numCols - 1);
+  for (unsigned i = 0; i < d; ++i)
     vertex.negateRow(i);
 
   return vertex;
@@ -229,7 +230,7 @@ mlir::presburger::detail::solveParametricEquations(FracMatrix equations) {
 /// we iterate over the vertex list, each time appending the vertex to the
 /// chambers where it is active and creating a new chamber if necessary.
 std::vector<std::pair<PresburgerRelation, std::vector<unsigned>>>
-mlir::presburger::detail::chamberDecomposition(
+mlir::presburger::detail::computeChamberDecomposition(
     std::vector<PresburgerRelation> activeRegions,
     std::vector<ParamPoint> vertices) {
   // We maintain a list of regions and their associated vertex sets,
@@ -428,7 +429,7 @@ mlir::presburger::detail::polytopeGeneratingFunction(PolyhedronH poly) {
   // they may share "faces" or "edges", but their intersection can only have
   // up to numVars-1 dimensions.
   std::vector<std::pair<PresburgerRelation, std::vector<unsigned>>> chambers =
-      chamberDecomposition(activeRegions, vertices);
+      computeChamberDecomposition(activeRegions, vertices);
 
   // Now, we compute the generating function. For each chamber, we iterate over
   // the vertices active in it, and compute the generating function for each
diff --git a/mlir/lib/Analysis/Presburger/IntegerRelation.cpp b/mlir/lib/Analysis/Presburger/IntegerRelation.cpp
index 5895b927435a424..201dc9bd4c8a1df 100644
--- a/mlir/lib/Analysis/Presburger/IntegerRelation.cpp
+++ b/mlir/lib/Analysis/Presburger/IntegerRelation.cpp
@@ -26,6 +26,7 @@
 #include "llvm/ADT/DenseMap.h"
 #include "llvm/ADT/DenseSet.h"
 #include "llvm/ADT/STLExtras.h"
+#include "llvm/ADT/Sequence.h"
 #include "llvm/ADT/SmallBitVector.h"
 #include "llvm/Support/Debug.h"
 #include "llvm/Support/raw_ostream.h"
@@ -2499,22 +2500,17 @@ void IntegerRelation::printSpace(raw_ostream &os) const {
 }
 
 void IntegerRelation::removeTrivialEqualities() {
-  std::vector<bool> isTrivialEquality;
-  for (int i = 0, e = getNumEqualities(); i < e; ++i) {
-    bool currentIsTrivial =
-        llvm::all_of(getEquality(i), [&](MPInt n) -> bool { return (n == 0); });
-    isTrivialEquality.push_back(currentIsTrivial);
-  }
-
-  unsigned pos = 0;
-  for (unsigned r = 0, e = getNumEqualities(); r < e; r++)
-    if (!isTrivialEquality[r])
-      equalities.copyRow(r, pos++);
-
-  equalities.resizeVertically(pos);
+  for (int i = getNumEqualities() - 1; i >= 0; --i)
+    if (rangeIsZero(getEquality(i)))
+      removeEquality(i);
 }
 
 bool IntegerRelation::isFullDim() {
+  if (getNumVars() == 0)
+    return true;
+  if (isEmpty())
+    return false;
+
   // If there is a non-trivial equality, the space cannot be full-dimensional.
   removeTrivialEqualities();
   if (getNumEqualities() > 0)
@@ -2523,24 +2519,11 @@ bool IntegerRelation::isFullDim() {
   // If along the direction of any of the inequalities, the upper and lower
   // optima are the same, then the region is not full-dimensional.
   Simplex simplex(*this);
-  for (unsigned i = 0, e = getNumInequalities(); i < e; i++) {
-    auto ineq = inequalities.getRow(i);
-    auto upOpt = simplex.computeOptimum(Simplex::Direction::Up, ineq);
-    auto downOpt = simplex.computeOptimum(Simplex::Direction::Down, ineq);
-
-    if (upOpt.getKind() == OptimumKind::Unbounded ||
-        downOpt.getKind() == OptimumKind::Unbounded)
-      continue;
-
-    // Check if the upper and lower optima are equal.
-    if (upOpt.getKind() == OptimumKind::Bounded &&
-        downOpt.getKind() == OptimumKind::Bounded && *upOpt == *downOpt)
-      return false;
-  }
-  // If none of the inequalities were such that the upper and lower optima
-  // along their direction were equal, then we conclude that the region is full
-  // dimensional.
-  return true;
+  bool fullDim = llvm::none_of(llvm::seq<int>(0, getNumInequalities()),
+                               [&](unsigned i) -> bool {
+                                 return simplex.isFlatAlong(getInequality(i));
+                               });
+  return fullDim;
 }
 
 void IntegerRelation::print(raw_ostream &os) const {
diff --git a/mlir/lib/Analysis/Presburger/Matrix.cpp b/mlir/lib/Analysis/Presburger/Matrix.cpp
index 35af156552c34d9..af4435bda4172bd 100644
--- a/mlir/lib/Analysis/Presburger/Matrix.cpp
+++ b/mlir/lib/Analysis/Presburger/Matrix.cpp
@@ -358,8 +358,8 @@ template <typename T>
 Matrix<T> Matrix<T>::getSubMatrix(unsigned fromRow, unsigned toRow,
                                   unsigned fromColumn,
                                   unsigned toColumn) const {
-  assert(toRow >= fromRow && "end of row range must be after beginning!");
-  assert(toColumn >= fromColumn && "end of row range must be after beginning!");
+  assert(fromRow <= toRow && "end of row range must be after beginning!");
+  assert(fromColumn <= toColumn && "end of row range must be after beginning!");
   Matrix<T> subMatrix(toRow - fromRow + 1, toColumn - fromColumn + 1);
   for (unsigned i = fromRow; i <= toRow; i++)
     for (unsigned j = fromColumn; j <= toColumn; j++)
diff --git a/mlir/lib/Analysis/Presburger/Simplex.cpp b/mlir/lib/Analysis/Presburger/Simplex.cpp
index f6d24b3f4dbd40e..ee3ec87e592c479 100644
--- a/mlir/lib/Analysis/Presburger/Simplex.cpp
+++ b/mlir/lib/Analysis/Presburger/Simplex.cpp
@@ -2104,12 +2104,20 @@ Simplex::computeIntegerBounds(ArrayRef<MPInt> coeffs) {
   return {minRoundedUp, maxRoundedDown};
 }
 
-bool Simplex::isValidEquality(ArrayRef<MPInt> coeffs) {
-  auto [downOpt, upOpt] = Simplex::computeIntegerBounds(coeffs);
-  if (upOpt.getKind() == OptimumKind::Bounded &&
-      downOpt.getKind() == OptimumKind::Bounded && *upOpt == *downOpt)
+bool Simplex::isFlatAlong(ArrayRef<MPInt> coeffs) {
+  assert(!isEmpty() && "cannot check for flatness of empty simplex!");
+  auto upOpt = computeOptimum(Simplex::Direction::Up, coeffs);
+  auto downOpt = computeOptimum(Simplex::Direction::Down, coeffs);
+
+  if (upOpt.getKind() != OptimumKind::Bounded)
     return false;
-  return true;
+  if (downOpt.getKind() != OptimumKind::Bounded)
+    return false;
+
+  // Check if the upper and lower optima are equal.
+  if (*upOpt == *downOpt)
+    return true;
+  return false;
 }
 
 void SimplexBase::print(raw_ostream &os) const {

>From 21086534cbb0c9fcc3e691be7983318300076937 Mon Sep 17 00:00:00 2001
From: Abhinav271828 <abhinav.m at research.iiit.ac.in>
Date: Wed, 24 Jan 2024 03:07:08 +0530
Subject: [PATCH 23/26] Minor

---
 mlir/lib/Analysis/Presburger/Matrix.cpp | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/mlir/lib/Analysis/Presburger/Matrix.cpp b/mlir/lib/Analysis/Presburger/Matrix.cpp
index af4435bda4172bd..0ba1dc410791403 100644
--- a/mlir/lib/Analysis/Presburger/Matrix.cpp
+++ b/mlir/lib/Analysis/Presburger/Matrix.cpp
@@ -361,8 +361,8 @@ Matrix<T> Matrix<T>::getSubMatrix(unsigned fromRow, unsigned toRow,
   assert(fromRow <= toRow && "end of row range must be after beginning!");
   assert(fromColumn <= toColumn && "end of row range must be after beginning!");
   Matrix<T> subMatrix(toRow - fromRow + 1, toColumn - fromColumn + 1);
-  for (unsigned i = fromRow; i <= toRow; i++)
-    for (unsigned j = fromColumn; j <= toColumn; j++)
+  for (unsigned i = fromRow; i <= toRow; ++i)
+    for (unsigned j = fromColumn; j <= toColumn; ++j)
       subMatrix(i - fromRow, j - fromColumn) = at(i, j);
   return subMatrix;
 }
@@ -388,7 +388,7 @@ Matrix<T>::splitByBitset(std::bitset<16> indicator) {
     else
       rowsForZero.appendExtraRow(getRow(i));
   }
-  return std::make_pair(rowsForOne, rowsForZero);
+  return {rowsForOne, rowsForZero};
 }
 
 template <typename T>

>From 540a0c15b660e7200c535ec15066a1cd25466bb1 Mon Sep 17 00:00:00 2001
From: Abhinav271828 <abhinav.m at research.iiit.ac.in>
Date: Wed, 24 Jan 2024 22:56:20 +0530
Subject: [PATCH 24/26] Minor fixes

---
 .../mlir/Analysis/Presburger/Barvinok.h       |  5 +-
 .../Analysis/Presburger/IntegerRelation.h     |  4 +-
 .../include/mlir/Analysis/Presburger/Matrix.h |  2 +-
 .../Analysis/Presburger/PresburgerRelation.h  |  4 ++
 .../mlir/Analysis/Presburger/Simplex.h        |  3 +-
 mlir/lib/Analysis/Presburger/Barvinok.cpp     | 51 ++++++++++---------
 .../Analysis/Presburger/IntegerRelation.cpp   |  8 ++-
 mlir/lib/Analysis/Presburger/Matrix.cpp       |  5 +-
 .../Presburger/PresburgerRelation.cpp         |  7 +++
 mlir/lib/Analysis/Presburger/Simplex.cpp      |  8 ++-
 10 files changed, 57 insertions(+), 40 deletions(-)

diff --git a/mlir/include/mlir/Analysis/Presburger/Barvinok.h b/mlir/include/mlir/Analysis/Presburger/Barvinok.h
index 4d3f643187247f2..f322842e9aff3c0 100644
--- a/mlir/include/mlir/Analysis/Presburger/Barvinok.h
+++ b/mlir/include/mlir/Analysis/Presburger/Barvinok.h
@@ -101,11 +101,12 @@ std::optional<ParamPoint> solveParametricEquations(FracMatrix equations);
 /// In the return type, the vertices are stored by their index in the
 /// `vertices` argument, i.e., the set {2, 3, 4} represents the vertex set
 /// {vertices[2], vertices[3], vertices[4]}.
+/// The ith relation corresponds to the activity region of the ith vertex set.
 /// Note that here, by disjoint, we mean that the intersection is not
 /// full-dimensional.
 std::vector<std::pair<PresburgerRelation, std::vector<unsigned>>>
-computeChamberDecomposition(std::vector<PresburgerRelation> activeRegions,
-                            std::vector<ParamPoint> vertices);
+computeChamberDecomposition(ArrayRef<PresburgerRelation> activeRegions,
+                            ArrayRef<ParamPoint> vertices);
 
 /// Compute the generating function corresponding to a polytope.
 /// All tangent cones of the polytope must be unimodular.
diff --git a/mlir/include/mlir/Analysis/Presburger/IntegerRelation.h b/mlir/include/mlir/Analysis/Presburger/IntegerRelation.h
index 3e76e4bd2cc0ce2..4c4576bf6d0d351 100644
--- a/mlir/include/mlir/Analysis/Presburger/IntegerRelation.h
+++ b/mlir/include/mlir/Analysis/Presburger/IntegerRelation.h
@@ -716,8 +716,10 @@ class IntegerRelation {
 
   // Verify whether the relation is full-dimensional, i.e.,
   // no equality holds for the relation.
+  //
   // If there are no variables, it always returns true.
-  // If there is at least one variable and the space is empty, it returns false.
+  // If there is at least one variable and the relation is empty, it returns
+  // false.
   bool isFullDim();
 
   void print(raw_ostream &os) const;
diff --git a/mlir/include/mlir/Analysis/Presburger/Matrix.h b/mlir/include/mlir/Analysis/Presburger/Matrix.h
index ae781a93a036ff6..9cd8cebeb2c84f4 100644
--- a/mlir/include/mlir/Analysis/Presburger/Matrix.h
+++ b/mlir/include/mlir/Analysis/Presburger/Matrix.h
@@ -311,7 +311,7 @@ class FracMatrix : public Matrix<Fraction> {
   // calls `y`, usually 3/4.
   void LLL(Fraction delta);
 
-  // Multiply each row of the matrix by the LCD of the denominators, thereby
+  // Multiply each row of the matrix by the LCM of the denominators, thereby
   // converting it to an integer matrix.
   IntMatrix normalizeRows();
 };
diff --git a/mlir/include/mlir/Analysis/Presburger/PresburgerRelation.h b/mlir/include/mlir/Analysis/Presburger/PresburgerRelation.h
index c6b00eca90733a0..35c89bc6e68466a 100644
--- a/mlir/include/mlir/Analysis/Presburger/PresburgerRelation.h
+++ b/mlir/include/mlir/Analysis/Presburger/PresburgerRelation.h
@@ -217,6 +217,10 @@ class PresburgerRelation {
   /// redundencies.
   PresburgerRelation simplify() const;
 
+  /// Return true if any of the disjuncts is full-dimensional
+  /// (see IntegerRelation::isFullDim()).
+  bool isFullDim() const;
+
   /// Print the set's internal state.
   void print(raw_ostream &os) const;
   void dump() const;
diff --git a/mlir/include/mlir/Analysis/Presburger/Simplex.h b/mlir/include/mlir/Analysis/Presburger/Simplex.h
index 8be9942a3dab922..7ee74c150867c1f 100644
--- a/mlir/include/mlir/Analysis/Presburger/Simplex.h
+++ b/mlir/include/mlir/Analysis/Presburger/Simplex.h
@@ -773,8 +773,9 @@ class Simplex : public SimplexBase {
 
   /// Check if the simplex takes only one rational value along the
   /// direction of `coeffs`.
+  ///
   /// `this` must be nonempty.
-  bool isFlatAlong(const ArrayRef<MPInt> coeffs);
+  bool isFlatAlong(ArrayRef<MPInt> coeffs);
 
   /// Returns true if the polytope is unbounded, i.e., extends to infinity in
   /// some direction. Otherwise, returns false.
diff --git a/mlir/lib/Analysis/Presburger/Barvinok.cpp b/mlir/lib/Analysis/Presburger/Barvinok.cpp
index 70ec9da4d530f3d..1d22880393846be 100644
--- a/mlir/lib/Analysis/Presburger/Barvinok.cpp
+++ b/mlir/lib/Analysis/Presburger/Barvinok.cpp
@@ -172,6 +172,8 @@ mlir::presburger::detail::solveParametricEquations(FracMatrix equations) {
           .determinant() == 0)
     return std::nullopt;
 
+  // Perform row operations to make each column all zeros except for the
+  // diagonal element, which is made to be one.
   for (unsigned i = 0; i < d; ++i) {
     // First ensure that the diagonal element is nonzero, by swapping
     // it with a row that is non-zero at column i.
@@ -196,7 +198,7 @@ mlir::presburger::detail::solveParametricEquations(FracMatrix equations) {
       // ith row, appropriately scaled.
       Fraction currentElement = equations(j, i);
       equations.addToRow(/*sourceRow=*/i, /*targetRow=*/j,
-                         -currentElement / diagElement);
+                         /*scale=*/-currentElement / diagElement);
     }
   }
 
@@ -229,10 +231,14 @@ mlir::presburger::detail::solveParametricEquations(FracMatrix equations) {
 /// We maintain a list of pairwise disjoint chambers and their vertex-sets;
 /// we iterate over the vertex list, each time appending the vertex to the
 /// chambers where it is active and creating a new chamber if necessary.
+///
+/// Given the region each vertex is active in, for each subset of vertices,
+/// the region that precisely this subset is in, is the intersection of the
+/// regions that these are active in, intersected with the complements of the
+/// remaining regions.
 std::vector<std::pair<PresburgerRelation, std::vector<unsigned>>>
 mlir::presburger::detail::computeChamberDecomposition(
-    std::vector<PresburgerRelation> activeRegions,
-    std::vector<ParamPoint> vertices) {
+    ArrayRef<PresburgerRelation> activeRegions, ArrayRef<ParamPoint> vertices) {
   // We maintain a list of regions and their associated vertex sets,
   // initialized with the first vertex and its corresponding activity region.
   std::vector<std::pair<PresburgerRelation, std::vector<unsigned>>> chambers = {
@@ -241,8 +247,6 @@ mlir::presburger::detail::computeChamberDecomposition(
   // of indices. Thus the set {2, 3, 4} represents the vertex set
   // {vertices[2], vertices[3], vertices[4]}.
 
-  std::vector<std::pair<PresburgerRelation, std::vector<unsigned>>> newChambers;
-
   // We iterate over the vertex set.
   // For each vertex v_j and its activity region R_j,
   // we examine all the current chambers R_i.
@@ -253,37 +257,35 @@ mlir::presburger::detail::computeChamberDecomposition(
   // Once we have examined all R_i, we add a final chamber
   // R_j - (union of all existing chambers),
   // in which only v_j is active.
+
+  // At each step, we define a new chamber list after considering vertex v_j,
+  // replacing and appending chambers as discussed above.
+  std::vector<std::pair<PresburgerRelation, std::vector<unsigned>>> newChambers;
   for (unsigned j = 1, e = vertices.size(); j < e; j++) {
     newChambers.clear();
 
     PresburgerRelation newRegion = activeRegions[j];
-    ParamPoint newVertex = vertices[j];
-
-    for (unsigned i = 0, f = chambers.size(); i < f; i++) {
-      auto [currentRegion, currentVertices] = chambers[i];
 
+    for (auto [currentRegion, currentVertices] : chambers) {
       // First, we check if the intersection of R_j and R_i.
       // It is a disjoint union of convex regions in the parameter space,
       // and so we know that it is full-dimensional if any of the disjuncts
       // is full-dimensional.
       PresburgerRelation intersection = currentRegion.intersect(newRegion);
-      bool isFullDim = intersection.getNumRangeVars() == 0 ||
-                       llvm::any_of(intersection.getAllDisjuncts(),
-                                    [&](IntegerRelation disjunct) -> bool {
-                                      return disjunct.isFullDim();
-                                    });
+      bool isFullDim =
+          intersection.getNumRangeVars() == 0 || intersection.isFullDim();
 
       // If the intersection is not full-dimensional, we do not modify
       // the chamber list.
       if (!isFullDim)
-        newChambers.push_back(chambers[i]);
+        newChambers.emplace_back(currentRegion, currentVertices);
       else {
         // If it is, we add the intersection and the difference as new chambers.
         PresburgerRelation subtraction = currentRegion.subtract(newRegion);
-        newChambers.push_back(std::make_pair(subtraction, currentVertices));
+        newChambers.emplace_back(subtraction, currentVertices);
 
         currentVertices.push_back(j);
-        newChambers.push_back(std::make_pair(intersection, currentVertices));
+        newChambers.emplace_back(intersection, currentVertices);
       }
     }
 
@@ -441,27 +443,28 @@ mlir::presburger::detail::polytopeGeneratingFunction(PolyhedronH poly) {
     auto [currentRegion, currentVertices] = chamber;
     GeneratingFunction chamberGf(numSymbols, {}, {}, {});
     for (unsigned i : currentVertices) {
-      // We collect the inequalities corresponding to each vertex.
+      // We collect the inequalities corresponding to each vertex to compute
+      // the tangent cone at that vertex.
       // We only need the coefficients of the variables (NOT the parameters)
       // as the generating function only depends on these.
       // We translate the cones to be pointed at the origin by making the
       // constant terms zero.
-      ConeH tgtCone = defineHRep(numVars);
-      for (unsigned j = 0; j < numVars; j++) {
-        for (unsigned k = 0; k < numVars; k++)
+      ConeH tangentCone = defineHRep(numVars);
+      for (unsigned j = 0; j < numVars; ++j) {
+        for (unsigned k = 0; k < numVars; ++k)
           ineq[k] = subsets[i](j, k);
-        tgtCone.addInequality(ineq);
+        tangentCone.addInequality(ineq);
       }
       // We assume that the tangent cone is unimodular.
       SmallVector<std::pair<int, ConeH>, 4> unimodCones = {
-          std::make_pair(1, tgtCone)};
+          std::make_pair(1, tangentCone)};
       for (std::pair<int, ConeH> signedCone : unimodCones) {
         auto [sign, cone] = signedCone;
         chamberGf = chamberGf +
                     unimodularConeGeneratingFunction(vertices[i], sign, cone);
       }
     }
-    gf.push_back(std::make_pair(currentRegion, chamberGf));
+    gf.emplace_back(currentRegion, chamberGf);
   }
   return gf;
 }
diff --git a/mlir/lib/Analysis/Presburger/IntegerRelation.cpp b/mlir/lib/Analysis/Presburger/IntegerRelation.cpp
index 201dc9bd4c8a1df..ab50f0dbec2cf3d 100644
--- a/mlir/lib/Analysis/Presburger/IntegerRelation.cpp
+++ b/mlir/lib/Analysis/Presburger/IntegerRelation.cpp
@@ -2519,11 +2519,9 @@ bool IntegerRelation::isFullDim() {
   // If along the direction of any of the inequalities, the upper and lower
   // optima are the same, then the region is not full-dimensional.
   Simplex simplex(*this);
-  bool fullDim = llvm::none_of(llvm::seq<int>(0, getNumInequalities()),
-                               [&](unsigned i) -> bool {
-                                 return simplex.isFlatAlong(getInequality(i));
-                               });
-  return fullDim;
+  return llvm::none_of(llvm::seq<int>(0, getNumInequalities()), [&](int i) {
+    return simplex.isFlatAlong(getInequality(i));
+  });
 }
 
 void IntegerRelation::print(raw_ostream &os) const {
diff --git a/mlir/lib/Analysis/Presburger/Matrix.cpp b/mlir/lib/Analysis/Presburger/Matrix.cpp
index 0ba1dc410791403..6780288463a7d8d 100644
--- a/mlir/lib/Analysis/Presburger/Matrix.cpp
+++ b/mlir/lib/Analysis/Presburger/Matrix.cpp
@@ -359,7 +359,10 @@ Matrix<T> Matrix<T>::getSubMatrix(unsigned fromRow, unsigned toRow,
                                   unsigned fromColumn,
                                   unsigned toColumn) const {
   assert(fromRow <= toRow && "end of row range must be after beginning!");
-  assert(fromColumn <= toColumn && "end of row range must be after beginning!");
+  assert(toRow < nRows && "end of row range out of bounds!");
+  assert(fromColumn <= toColumn &&
+         "end of column range must be after beginning!");
+  assert(toColumn < nColumns && "end of column range out of bounds!");
   Matrix<T> subMatrix(toRow - fromRow + 1, toColumn - fromColumn + 1);
   for (unsigned i = fromRow; i <= toRow; ++i)
     for (unsigned j = fromColumn; j <= toColumn; ++j)
diff --git a/mlir/lib/Analysis/Presburger/PresburgerRelation.cpp b/mlir/lib/Analysis/Presburger/PresburgerRelation.cpp
index 787fc1c659a12e3..ee56ce696e86e06 100644
--- a/mlir/lib/Analysis/Presburger/PresburgerRelation.cpp
+++ b/mlir/lib/Analysis/Presburger/PresburgerRelation.cpp
@@ -1041,6 +1041,13 @@ PresburgerRelation PresburgerRelation::simplify() const {
   return result;
 }
 
+bool PresburgerRelation::isFullDim() const {
+  bool fullDim = llvm::any_of(getAllDisjuncts(), [&](IntegerRelation disjunct) {
+    return disjunct.isFullDim();
+  });
+  return fullDim;
+}
+
 void PresburgerRelation::print(raw_ostream &os) const {
   os << "Number of Disjuncts: " << getNumDisjuncts() << "\n";
   for (const IntegerRelation &disjunct : disjuncts) {
diff --git a/mlir/lib/Analysis/Presburger/Simplex.cpp b/mlir/lib/Analysis/Presburger/Simplex.cpp
index ee3ec87e592c479..592d96ee766fb06 100644
--- a/mlir/lib/Analysis/Presburger/Simplex.cpp
+++ b/mlir/lib/Analysis/Presburger/Simplex.cpp
@@ -2109,15 +2109,13 @@ bool Simplex::isFlatAlong(ArrayRef<MPInt> coeffs) {
   auto upOpt = computeOptimum(Simplex::Direction::Up, coeffs);
   auto downOpt = computeOptimum(Simplex::Direction::Down, coeffs);
 
-  if (upOpt.getKind() != OptimumKind::Bounded)
+  if (!upOpt.isBounded())
     return false;
-  if (downOpt.getKind() != OptimumKind::Bounded)
+  if (!downOpt.isBounded())
     return false;
 
   // Check if the upper and lower optima are equal.
-  if (*upOpt == *downOpt)
-    return true;
-  return false;
+  return *upOpt == *downOpt;
 }
 
 void SimplexBase::print(raw_ostream &os) const {

>From 5dc100f478a14635a41f11ac827ed7212651ca38 Mon Sep 17 00:00:00 2001
From: Abhinav271828 <abhinav.m at research.iiit.ac.in>
Date: Wed, 24 Jan 2024 22:58:26 +0530
Subject: [PATCH 25/26] Use const

---
 mlir/include/mlir/Analysis/Presburger/Matrix.h | 2 +-
 mlir/lib/Analysis/Presburger/Matrix.cpp        | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/mlir/include/mlir/Analysis/Presburger/Matrix.h b/mlir/include/mlir/Analysis/Presburger/Matrix.h
index 9cd8cebeb2c84f4..fe260568a22d874 100644
--- a/mlir/include/mlir/Analysis/Presburger/Matrix.h
+++ b/mlir/include/mlir/Analysis/Presburger/Matrix.h
@@ -313,7 +313,7 @@ class FracMatrix : public Matrix<Fraction> {
 
   // Multiply each row of the matrix by the LCM of the denominators, thereby
   // converting it to an integer matrix.
-  IntMatrix normalizeRows();
+  IntMatrix normalizeRows() const;
 };
 
 } // namespace presburger
diff --git a/mlir/lib/Analysis/Presburger/Matrix.cpp b/mlir/lib/Analysis/Presburger/Matrix.cpp
index 6780288463a7d8d..6d734c232f18281 100644
--- a/mlir/lib/Analysis/Presburger/Matrix.cpp
+++ b/mlir/lib/Analysis/Presburger/Matrix.cpp
@@ -729,7 +729,7 @@ void FracMatrix::LLL(Fraction delta) {
   }
 }
 
-IntMatrix FracMatrix::normalizeRows() {
+IntMatrix FracMatrix::normalizeRows() const {
   unsigned numRows = getNumRows();
   unsigned numColumns = getNumColumns();
   IntMatrix normalized(numRows, numColumns);

>From aaedcc0c149cce62917f5d9eaa1626041f57fea3 Mon Sep 17 00:00:00 2001
From: Abhinav271828 <abhinav.m at research.iiit.ac.in>
Date: Wed, 24 Jan 2024 23:01:12 +0530
Subject: [PATCH 26/26] Rename

---
 mlir/include/mlir/Analysis/Presburger/Barvinok.h    | 7 ++++---
 mlir/lib/Analysis/Presburger/Barvinok.cpp           | 9 +++++----
 mlir/unittests/Analysis/Presburger/BarvinokTest.cpp | 9 +++++----
 3 files changed, 14 insertions(+), 11 deletions(-)

diff --git a/mlir/include/mlir/Analysis/Presburger/Barvinok.h b/mlir/include/mlir/Analysis/Presburger/Barvinok.h
index f322842e9aff3c0..729a49b0e8e3294 100644
--- a/mlir/include/mlir/Analysis/Presburger/Barvinok.h
+++ b/mlir/include/mlir/Analysis/Presburger/Barvinok.h
@@ -85,8 +85,9 @@ ConeH getDual(ConeV cone);
 
 /// Compute the generating function for a unimodular cone.
 /// The input cone must be unimodular; it assert-fails otherwise.
-GeneratingFunction unimodularConeGeneratingFunction(ParamPoint vertex, int sign,
-                                                    ConeH cone);
+GeneratingFunction computeUnimodularConeGeneratingFunction(ParamPoint vertex,
+                                                           int sign,
+                                                           ConeH cone);
 
 /// Find the solution of a set of equations that express affine constraints
 /// between a set of variables and a set of parameters. The solution expresses
@@ -111,7 +112,7 @@ computeChamberDecomposition(ArrayRef<PresburgerRelation> activeRegions,
 /// Compute the generating function corresponding to a polytope.
 /// All tangent cones of the polytope must be unimodular.
 std::vector<std::pair<PresburgerRelation, GeneratingFunction>>
-polytopeGeneratingFunction(PolyhedronH poly);
+computePolytopeGeneratingFunction(PolyhedronH poly);
 
 /// Find a vector that is not orthogonal to any of the given vectors,
 /// i.e., has nonzero dot product with those of the given vectors
diff --git a/mlir/lib/Analysis/Presburger/Barvinok.cpp b/mlir/lib/Analysis/Presburger/Barvinok.cpp
index 1d22880393846be..b997b48d960e38a 100644
--- a/mlir/lib/Analysis/Presburger/Barvinok.cpp
+++ b/mlir/lib/Analysis/Presburger/Barvinok.cpp
@@ -77,7 +77,8 @@ MPInt mlir::presburger::detail::getIndex(ConeV cone) {
 /// num is computed by expressing the vertex as a weighted
 /// sum of the generators, and then taking the floor of the
 /// coefficients.
-GeneratingFunction mlir::presburger::detail::unimodularConeGeneratingFunction(
+GeneratingFunction
+mlir::presburger::detail::computeUnimodularConeGeneratingFunction(
     ParamPoint vertex, int sign, ConeH cone) {
   // Consider a cone with H-representation [0  -1].
   //                                       [-1 -2]
@@ -328,7 +329,7 @@ mlir::presburger::detail::computeChamberDecomposition(
 /// polytopes using Barvinok's rational functions." Algorithmica 48 (2007):
 /// 37-66.
 std::vector<std::pair<PresburgerRelation, GeneratingFunction>>
-mlir::presburger::detail::polytopeGeneratingFunction(PolyhedronH poly) {
+mlir::presburger::detail::computePolytopeGeneratingFunction(PolyhedronH poly) {
   unsigned numVars = poly.getNumRangeVars();
   unsigned numSymbols = poly.getNumSymbolVars();
   unsigned numIneqs = poly.getNumInequalities();
@@ -460,8 +461,8 @@ mlir::presburger::detail::polytopeGeneratingFunction(PolyhedronH poly) {
           std::make_pair(1, tangentCone)};
       for (std::pair<int, ConeH> signedCone : unimodCones) {
         auto [sign, cone] = signedCone;
-        chamberGf = chamberGf +
-                    unimodularConeGeneratingFunction(vertices[i], sign, cone);
+        chamberGf = chamberGf + computeUnimodularConeGeneratingFunction(
+                                    vertices[i], sign, cone);
       }
     }
     gf.emplace_back(currentRegion, chamberGf);
diff --git a/mlir/unittests/Analysis/Presburger/BarvinokTest.cpp b/mlir/unittests/Analysis/Presburger/BarvinokTest.cpp
index 09301d5e54e0e0c..50a6cfa0fb6ccb4 100644
--- a/mlir/unittests/Analysis/Presburger/BarvinokTest.cpp
+++ b/mlir/unittests/Analysis/Presburger/BarvinokTest.cpp
@@ -59,7 +59,8 @@ TEST(BarvinokTest, unimodularConeGeneratingFunction) {
   ParamPoint vertex =
       makeFracMatrix(2, 3, {{2, 2, 0}, {-1, -Fraction(1, 2), 1}});
 
-  GeneratingFunction gf = unimodularConeGeneratingFunction(vertex, 1, cone);
+  GeneratingFunction gf =
+      computeUnimodularConeGeneratingFunction(vertex, 1, cone);
 
   EXPECT_EQ_REPR_GENERATINGFUNCTION(
       gf, GeneratingFunction(
@@ -74,7 +75,7 @@ TEST(BarvinokTest, unimodularConeGeneratingFunction) {
 
   vertex = makeFracMatrix(3, 2, {{5, 2}, {6, 2}, {7, 1}});
 
-  gf = unimodularConeGeneratingFunction(vertex, 1, cone);
+  gf = computeUnimodularConeGeneratingFunction(vertex, 1, cone);
 
   EXPECT_EQ_REPR_GENERATINGFUNCTION(
       gf,
@@ -252,7 +253,7 @@ TEST(BarvinokTest, computeNumTermsPolytope) {
     poly.addInequality(ineqs.getRow(i));
 
   std::vector<std::pair<PresburgerRelation, GeneratingFunction>> count =
-      polytopeGeneratingFunction(poly);
+      computePolytopeGeneratingFunction(poly);
   // There is only one chamber, as it is non-parametric.
   EXPECT_EQ(count.size(), 1u);
 
@@ -281,7 +282,7 @@ TEST(BarvinokTest, computeNumTermsPolytope) {
   for (unsigned i = 0; i < 3; i++)
     poly.addInequality(ineqs.getRow(i));
 
-  count = polytopeGeneratingFunction(poly);
+  count = computePolytopeGeneratingFunction(poly);
   // There is only one chamber: p ≥ 0
   EXPECT_EQ(count.size(), 1u);
 



More information about the Mlir-commits mailing list