[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
Tue Jan 23 05:15:23 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/15] 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/15] 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/15] 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/15] 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/15] 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/15] 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/15] ../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/15] 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/15] 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/15] 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/15] 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/15] 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/15] 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/15] 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/15] 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);
+
+}



More information about the Mlir-commits mailing list