[Mlir-commits] [mlir] [MLIR][Presburger] Implement vertex enumeration and chamber decomposition for polytope generating function computation. (PR #78987)
llvmlistbot at llvm.org
llvmlistbot at llvm.org
Tue Jan 30 09:49: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/41] 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 b70ec33b69323..6feb7e8f09233 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 c38eab6efd0fc..35f15319b4568 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 c476a022a4827..e7b03c44d33e3 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 e0fd0dd8caa4d..c8cec45a21dc4 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 7d2a63d17676f..578b95571fff2 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/41] 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 6feb7e8f09233..1fafa1e709d5e 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 c8cec45a21dc4..f9f75a8a0dcef 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 578b95571fff2..edb95b4072e1d 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/41] 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 1fafa1e709d5e..320848348552c 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 f9f75a8a0dcef..3753f715c7d92 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/41] 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 320848348552c..9c940231d5e14 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 3753f715c7d92..ff202984bd336 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/41] 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 9c940231d5e14..a5bc49d79408c 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 ff202984bd336..3afc9ad8638cc 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/41] 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 e7b03c44d33e3..141dc1f8fa97e 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 edb95b4072e1d..b05075b2892ac 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/41] ../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 b05075b2892ac..1683f5cf41aee 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/41] 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 3afc9ad8638cc..1a3dd9a32e50a 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/41] 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 1a3dd9a32e50a..b3e68d1f65f25 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/41] 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 a5bc49d79408c..b4d1d6de53079 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 141dc1f8fa97e..b60d68aab4cf4 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 9482f69b31cd6..1939a41e52bc6 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 b3e68d1f65f25..6452b2d37838d 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 42bbc3363d583..f6d24b3f4dbd4 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/41] 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 0d4a593a95b1c..4a5bba17e2045 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 6452b2d37838d..67af0eb452bf3 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 bd7f7f58a932f..be8dac636e0af 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/41] 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 67af0eb452bf3..8c4de595b685c 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/41] 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 b4d1d6de53079..34055e93242c4 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 4a5bba17e2045..ae781a93a036f 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 8c4de595b685c..00a1ee97cd38d 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 be8dac636e0af..099cb6fa09e76 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/41] 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 34055e93242c4..558498e991040 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 00a1ee97cd38d..4b32489943bb0 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 099cb6fa09e76..35af156552c34 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/41] 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 4b32489943bb0..82fbee5a05b44 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 919aaa7a42859..9770d03c90a9f 100644
--- a/mlir/unittests/Analysis/Presburger/BarvinokTest.cpp
+++ b/mlir/unittests/Analysis/Presburger/BarvinokTest.cpp
@@ -125,7 +125,7 @@ TEST(BarvinokTest, getCoefficientInRationalFunction) {
EXPECT_EQ(coeff.getConstantTerm(), Fraction(55, 64));
}
-TEST(BarvinokTest, computeNumTerms) {
+TEST(BarvinokTest, computeNumTermsCone) {
// The following test is taken from
// Verdoolaege, Sven, et al. "Counting integer points in parametric
// polytopes using Barvinok's rational functions." Algorithmica 48 (2007):
@@ -233,4 +233,20 @@ TEST(BarvinokTest, computeNumTerms) {
for (unsigned j = 0; j < 2; j++)
for (unsigned k = 0; k < 2; k++)
EXPECT_EQ(count[i][j][k], 1);
-}
\ No newline at end of file
+}
+
+TEST(BarvinokTest, computeNumTermsPolytope) {
+ IntMatrix ineqs = makeIntMatrix(6, 4, {{1, 0, 0, 0},
+ {0, 1, 0, 0},
+ {0, 0, 1, 0},
+ {-1, 0, 0, 1},
+ {0, -1, 0, 1},
+ {0, 0, -1, 1}});
+ PolyhedronH poly = defineHRep(3);
+ for (unsigned i = 0; i < 6; i++)
+ poly.addInequality(ineqs.getRow(i));
+
+ std::vector<std::pair<PresburgerRelation, GeneratingFunction>> count = polytopeGeneratingFunction(poly);
+ EXPECT_EQ(count.size(), 1u);
+
+}
>From fe73f449fd801c140902481515bd1e6a0c825fca Mon Sep 17 00:00:00 2001
From: Abhinav271828 <abhinav.m at research.iiit.ac.in>
Date: Tue, 23 Jan 2024 18:48:49 +0530
Subject: [PATCH 16/41] Minor fix
---
mlir/lib/Analysis/Presburger/Barvinok.cpp | 5 +++--
1 file changed, 3 insertions(+), 2 deletions(-)
diff --git a/mlir/lib/Analysis/Presburger/Barvinok.cpp b/mlir/lib/Analysis/Presburger/Barvinok.cpp
index 82fbee5a05b44..ae66bedd21cd2 100644
--- a/mlir/lib/Analysis/Presburger/Barvinok.cpp
+++ b/mlir/lib/Analysis/Presburger/Barvinok.cpp
@@ -85,7 +85,7 @@ GeneratingFunction mlir::presburger::detail::unimodularConeGeneratingFunction(
// [-1 -1/2 1]
// `cone` must be unimodular.
- assert(getIndex(getDual(cone)) == 1 && "input cone is not unimodular!");
+ assert(abs(getIndex(getDual(cone))) == 1 && "input cone is not unimodular!");
unsigned numVar = cone.getNumVars();
unsigned numIneq = cone.getNumInequalities();
@@ -443,11 +443,12 @@ mlir::presburger::detail::polytopeGeneratingFunction(PolyhedronH poly) {
// We collect the inequalities corresponding to each vertex.
// We only need the coefficients of the variables (NOT the parameters)
// as the generating function only depends on these.
+ // We translate the cones to be pointed at the origin by making the
+ // constant terms zero.
ConeH tgtCone = defineHRep(numVars);
for (unsigned j = 0; j < numVars; j++) {
for (unsigned k = 0; k < numVars; k++)
ineq[k] = subsets[i](j, k);
- ineq[numVars] = subsets[i](j, numVars + numSymbols);
tgtCone.addInequality(ineq);
}
// We assume that the tangent cone is unimodular.
>From 2d8512f8cf5bb576f33380289051a97afd544b5c Mon Sep 17 00:00:00 2001
From: Abhinav271828 <abhinav.m at research.iiit.ac.in>
Date: Tue, 23 Jan 2024 19:32:06 +0530
Subject: [PATCH 17/41] Add test for pgf
---
.../Analysis/Presburger/BarvinokTest.cpp | 48 ++++++++++++++-----
1 file changed, 37 insertions(+), 11 deletions(-)
diff --git a/mlir/unittests/Analysis/Presburger/BarvinokTest.cpp b/mlir/unittests/Analysis/Presburger/BarvinokTest.cpp
index 9770d03c90a9f..cd616fb7bbf3c 100644
--- a/mlir/unittests/Analysis/Presburger/BarvinokTest.cpp
+++ b/mlir/unittests/Analysis/Presburger/BarvinokTest.cpp
@@ -235,18 +235,44 @@ TEST(BarvinokTest, computeNumTermsCone) {
EXPECT_EQ(count[i][j][k], 1);
}
+/// We define some simple polyhedra with unimodular tangent cones and verify
+/// that the returned generating functions correspond to those calculated by
+/// hand.
TEST(BarvinokTest, computeNumTermsPolytope) {
- IntMatrix ineqs = makeIntMatrix(6, 4, {{1, 0, 0, 0},
- {0, 1, 0, 0},
- {0, 0, 1, 0},
- {-1, 0, 0, 1},
- {0, -1, 0, 1},
- {0, 0, -1, 1}});
- PolyhedronH poly = defineHRep(3);
- for (unsigned i = 0; i < 6; i++)
- poly.addInequality(ineqs.getRow(i));
-
- std::vector<std::pair<PresburgerRelation, GeneratingFunction>> count = polytopeGeneratingFunction(poly);
+ IntMatrix ineqs = makeIntMatrix(6, 4,
+ {{1, 0, 0, 0},
+ {0, 1, 0, 0},
+ {0, 0, 1, 0},
+ {-1, 0, 0, 1},
+ {0, -1, 0, 1},
+ {0, 0, -1, 1}});
+ PolyhedronH poly = defineHRep(3);
+ for (unsigned i = 0; i < 6; i++)
+ poly.addInequality(ineqs.getRow(i));
+
+ std::vector<std::pair<PresburgerRelation, GeneratingFunction>> count =
+ polytopeGeneratingFunction(poly);
+ // There is only one chamber, as it is non-parametric.
EXPECT_EQ(count.size(), 1u);
+ GeneratingFunction gf = count[0].second;
+ EXPECT_EQ_REPR_GENERATINGFUNCTION(
+ gf,
+ GeneratingFunction(
+ 0, {1, 1, 1, 1, 1, 1, 1, 1},
+ {makeFracMatrix(1, 3, {{1, 1, 1}}), makeFracMatrix(1, 3, {{0, 1, 1}}),
+ makeFracMatrix(1, 3, {{0, 1, 1}}), makeFracMatrix(1, 3, {{0, 0, 1}}),
+ makeFracMatrix(1, 3, {{0, 1, 1}}), makeFracMatrix(1, 3, {{0, 0, 1}}),
+ makeFracMatrix(1, 3, {{0, 0, 1}}),
+ makeFracMatrix(1, 3, {{0, 0, 0}})},
+ {{{-1, 0, 0}, {0, -1, 0}, {0, 0, -1}},
+ {{1, 0, 0}, {0, -1, 0}, {0, 0, -1}},
+ {{0, 1, 0}, {-1, 0, 0}, {0, 0, -1}},
+ {{1, 0, 0}, {0, 1, 0}, {0, 0, -1}},
+ {{0, 0, 1}, {-1, 0, 0}, {0, -1, 0}},
+ {{1, 0, 0}, {0, 0, 1}, {0, -1, 0}},
+ {{0, 1, 0}, {0, 0, 1}, {-1, 0, 0}},
+ {{1, 0, 0}, {0, 1, 0}, {0, 0, 1}}}));
+
+
}
>From 7c2906dcb6150f78a89188d54bf5c9c7d7d7a9d2 Mon Sep 17 00:00:00 2001
From: Abhinav271828 <abhinav.m at research.iiit.ac.in>
Date: Tue, 23 Jan 2024 19:34:58 +0530
Subject: [PATCH 18/41] Fix doc
---
mlir/include/mlir/Analysis/Presburger/Barvinok.h | 10 ++++++----
1 file changed, 6 insertions(+), 4 deletions(-)
diff --git a/mlir/include/mlir/Analysis/Presburger/Barvinok.h b/mlir/include/mlir/Analysis/Presburger/Barvinok.h
index 558498e991040..e6d4af510b1ca 100644
--- a/mlir/include/mlir/Analysis/Presburger/Barvinok.h
+++ b/mlir/include/mlir/Analysis/Presburger/Barvinok.h
@@ -49,14 +49,16 @@ using PolyhedronV = IntMatrix;
using ConeH = PolyhedronH;
using ConeV = PolyhedronV;
-inline ConeH defineHRep(int numVars, int numSymbols = 0) {
+inline PolyhedronH defineHRep(int numVars, int numSymbols = 0) {
// We don't distinguish between domain and range variables, so
// we set the number of domain variables as 0 and the number of
// range variables as the number of actual variables.
- // There are no symbols (we don't work with parametric cones) and no local
- // (existentially quantified) variables.
+ // numSymbols is the number of parameters; it is
+ // There are no local (existentially quantified) variables.
+ // The number of symbols is the number of parameters. By default, we consider
+ // nonparametric polyhedra.
// Once the cone is defined, we use `addInequality()` to set inequalities.
- return ConeH(PresburgerSpace::getSetSpace(/*numDims=*/numVars,
+ return PolyhedronH(PresburgerSpace::getSetSpace(/*numDims=*/numVars,
/*numSymbols=*/numSymbols,
/*numLocals=*/0));
}
>From 8b1169659263ca186f8dfc2b68ad5ef9f343caad Mon Sep 17 00:00:00 2001
From: Abhinav271828 <abhinav.m at research.iiit.ac.in>
Date: Tue, 23 Jan 2024 19:58:04 +0530
Subject: [PATCH 19/41] Fix removeTrivialEquality()
---
.../Analysis/Presburger/IntegerRelation.cpp | 29 ++++++++++---------
1 file changed, 15 insertions(+), 14 deletions(-)
diff --git a/mlir/lib/Analysis/Presburger/IntegerRelation.cpp b/mlir/lib/Analysis/Presburger/IntegerRelation.cpp
index 1683f5cf41aee..5895b927435a4 100644
--- a/mlir/lib/Analysis/Presburger/IntegerRelation.cpp
+++ b/mlir/lib/Analysis/Presburger/IntegerRelation.cpp
@@ -2499,30 +2499,31 @@ void IntegerRelation::printSpace(raw_ostream &os) const {
}
void IntegerRelation::removeTrivialEqualities() {
- bool flag;
- for (unsigned i = 0, e = getNumInequalities(); i < e; i++) {
- flag = true;
- for (unsigned j = 0, f = getNumVars(); j < f + 1; j++)
- if (atEq(i, j) != 0)
- flag = false;
- if (flag)
- removeEquality(i);
+ std::vector<bool> isTrivialEquality;
+ for (int i = 0, e = getNumEqualities(); i < e; ++i) {
+ bool currentIsTrivial =
+ llvm::all_of(getEquality(i), [&](MPInt n) -> bool { return (n == 0); });
+ isTrivialEquality.push_back(currentIsTrivial);
}
-}
-bool IntegerRelation::isFullDim() {
- removeTrivialEqualities();
+ unsigned pos = 0;
+ for (unsigned r = 0, e = getNumEqualities(); r < e; r++)
+ if (!isTrivialEquality[r])
+ equalities.copyRow(r, pos++);
- unsigned e = getNumInequalities();
+ equalities.resizeVertically(pos);
+}
+bool IntegerRelation::isFullDim() {
// If there is a non-trivial equality, the space cannot be full-dimensional.
- if (e > 0)
+ removeTrivialEqualities();
+ if (getNumEqualities() > 0)
return false;
// If along the direction of any of the inequalities, the upper and lower
// optima are the same, then the region is not full-dimensional.
Simplex simplex(*this);
- for (unsigned i = 0; i < e; i++) {
+ for (unsigned i = 0, e = getNumInequalities(); i < e; i++) {
auto ineq = inequalities.getRow(i);
auto upOpt = simplex.computeOptimum(Simplex::Direction::Up, ineq);
auto downOpt = simplex.computeOptimum(Simplex::Direction::Down, ineq);
>From c9448e38f24f92a3752d01ce7715f1ff87229fe0 Mon Sep 17 00:00:00 2001
From: Abhinav271828 <abhinav.m at research.iiit.ac.in>
Date: Tue, 23 Jan 2024 20:03:52 +0530
Subject: [PATCH 20/41] Formatting
---
mlir/include/mlir/Analysis/Presburger/Barvinok.h | 4 ++--
mlir/unittests/Analysis/Presburger/BarvinokTest.cpp | 10 ++++++++++
2 files changed, 12 insertions(+), 2 deletions(-)
diff --git a/mlir/include/mlir/Analysis/Presburger/Barvinok.h b/mlir/include/mlir/Analysis/Presburger/Barvinok.h
index e6d4af510b1ca..3a80d201730f5 100644
--- a/mlir/include/mlir/Analysis/Presburger/Barvinok.h
+++ b/mlir/include/mlir/Analysis/Presburger/Barvinok.h
@@ -59,8 +59,8 @@ inline PolyhedronH defineHRep(int numVars, int numSymbols = 0) {
// nonparametric polyhedra.
// Once the cone is defined, we use `addInequality()` to set inequalities.
return PolyhedronH(PresburgerSpace::getSetSpace(/*numDims=*/numVars,
- /*numSymbols=*/numSymbols,
- /*numLocals=*/0));
+ /*numSymbols=*/numSymbols,
+ /*numLocals=*/0));
}
/// Get the index of a cone, i.e., the volume of the parallelepiped
diff --git a/mlir/unittests/Analysis/Presburger/BarvinokTest.cpp b/mlir/unittests/Analysis/Presburger/BarvinokTest.cpp
index cd616fb7bbf3c..81c9b79bb14a3 100644
--- a/mlir/unittests/Analysis/Presburger/BarvinokTest.cpp
+++ b/mlir/unittests/Analysis/Presburger/BarvinokTest.cpp
@@ -239,6 +239,7 @@ TEST(BarvinokTest, computeNumTermsCone) {
/// that the returned generating functions correspond to those calculated by
/// hand.
TEST(BarvinokTest, computeNumTermsPolytope) {
+ // A cube of side 1.
IntMatrix ineqs = makeIntMatrix(6, 4,
{{1, 0, 0, 0},
{0, 1, 0, 0},
@@ -274,5 +275,14 @@ TEST(BarvinokTest, computeNumTermsPolytope) {
{{0, 1, 0}, {0, 0, 1}, {-1, 0, 0}},
{{1, 0, 0}, {0, 1, 0}, {0, 0, 1}}}));
+ // A right-angled triangle with side p.
+ ineqs = makeIntMatrix(3, 4, {{1, 0, 0, 0}, {0, 1, 0, 0}, {-1, -1, 1, 0}});
+ poly = defineHRep(2, 1);
+ for (unsigned i = 0; i < 3; i++)
+ poly.addInequality(ineqs.getRow(i));
+
+ count = polytopeGeneratingFunction(poly);
+ // There is only one chamber: p ≥ 0
+ EXPECT_EQ(count.size(), 1u);
}
>From 8a21ff3937ff0597514047318d0585e1225148de Mon Sep 17 00:00:00 2001
From: Abhinav271828 <abhinav.m at research.iiit.ac.in>
Date: Tue, 23 Jan 2024 20:25:30 +0530
Subject: [PATCH 21/41] Add test for pgf
---
mlir/unittests/Analysis/Presburger/BarvinokTest.cpp | 10 +++++++++-
1 file changed, 9 insertions(+), 1 deletion(-)
diff --git a/mlir/unittests/Analysis/Presburger/BarvinokTest.cpp b/mlir/unittests/Analysis/Presburger/BarvinokTest.cpp
index 81c9b79bb14a3..09301d5e54e0e 100644
--- a/mlir/unittests/Analysis/Presburger/BarvinokTest.cpp
+++ b/mlir/unittests/Analysis/Presburger/BarvinokTest.cpp
@@ -280,9 +280,17 @@ TEST(BarvinokTest, computeNumTermsPolytope) {
poly = defineHRep(2, 1);
for (unsigned i = 0; i < 3; i++)
poly.addInequality(ineqs.getRow(i));
-
+
count = polytopeGeneratingFunction(poly);
// There is only one chamber: p ≥ 0
EXPECT_EQ(count.size(), 1u);
+ gf = count[0].second;
+ EXPECT_EQ_REPR_GENERATINGFUNCTION(
+ gf, GeneratingFunction(
+ 1, {1, 1, 1},
+ {makeFracMatrix(2, 2, {{0, 1}, {0, 0}}),
+ makeFracMatrix(2, 2, {{0, 1}, {0, 0}}),
+ makeFracMatrix(2, 2, {{0, 0}, {0, 0}})},
+ {{{-1, 1}, {-1, 0}}, {{1, -1}, {0, -1}}, {{1, 0}, {0, 1}}}));
}
>From 57f4da54904e51850dff11e416abd13096b4baa6 Mon Sep 17 00:00:00 2001
From: Abhinav271828 <abhinav.m at research.iiit.ac.in>
Date: Wed, 24 Jan 2024 03:04:16 +0530
Subject: [PATCH 22/41] Minor fixes
---
.../mlir/Analysis/Presburger/Barvinok.h | 4 +-
.../Analysis/Presburger/GeneratingFunction.h | 2 +-
.../Analysis/Presburger/IntegerRelation.h | 2 +
.../mlir/Analysis/Presburger/Simplex.h | 7 +--
mlir/lib/Analysis/Presburger/Barvinok.cpp | 39 ++++++++--------
.../Analysis/Presburger/IntegerRelation.cpp | 45 ++++++-------------
mlir/lib/Analysis/Presburger/Matrix.cpp | 4 +-
mlir/lib/Analysis/Presburger/Simplex.cpp | 18 +++++---
8 files changed, 58 insertions(+), 63 deletions(-)
diff --git a/mlir/include/mlir/Analysis/Presburger/Barvinok.h b/mlir/include/mlir/Analysis/Presburger/Barvinok.h
index 3a80d201730f5..4d3f643187247 100644
--- a/mlir/include/mlir/Analysis/Presburger/Barvinok.h
+++ b/mlir/include/mlir/Analysis/Presburger/Barvinok.h
@@ -104,8 +104,8 @@ std::optional<ParamPoint> solveParametricEquations(FracMatrix equations);
/// Note that here, by disjoint, we mean that the intersection is not
/// full-dimensional.
std::vector<std::pair<PresburgerRelation, std::vector<unsigned>>>
-chamberDecomposition(std::vector<PresburgerRelation> activeRegions,
- std::vector<ParamPoint> vertices);
+computeChamberDecomposition(std::vector<PresburgerRelation> activeRegions,
+ std::vector<ParamPoint> vertices);
/// Compute the generating function corresponding to a polytope.
/// All tangent cones of the polytope must be unimodular.
diff --git a/mlir/include/mlir/Analysis/Presburger/GeneratingFunction.h b/mlir/include/mlir/Analysis/Presburger/GeneratingFunction.h
index 35f15319b4568..db5b6b6a95918 100644
--- a/mlir/include/mlir/Analysis/Presburger/GeneratingFunction.h
+++ b/mlir/include/mlir/Analysis/Presburger/GeneratingFunction.h
@@ -72,7 +72,7 @@ class GeneratingFunction {
return denominators;
}
- GeneratingFunction operator+(GeneratingFunction gf) const {
+ GeneratingFunction operator+(const GeneratingFunction &gf) const {
assert(numParam == gf.getNumParams() &&
"two generating functions with different numbers of parameters "
"cannot be added!");
diff --git a/mlir/include/mlir/Analysis/Presburger/IntegerRelation.h b/mlir/include/mlir/Analysis/Presburger/IntegerRelation.h
index b60d68aab4cf4..3e76e4bd2cc0c 100644
--- a/mlir/include/mlir/Analysis/Presburger/IntegerRelation.h
+++ b/mlir/include/mlir/Analysis/Presburger/IntegerRelation.h
@@ -716,6 +716,8 @@ class IntegerRelation {
// Verify whether the relation is full-dimensional, i.e.,
// no equality holds for the relation.
+ // If there are no variables, it always returns true.
+ // If there is at least one variable and the space is empty, it returns false.
bool isFullDim();
void print(raw_ostream &os) const;
diff --git a/mlir/include/mlir/Analysis/Presburger/Simplex.h b/mlir/include/mlir/Analysis/Presburger/Simplex.h
index 1939a41e52bc6..8be9942a3dab9 100644
--- a/mlir/include/mlir/Analysis/Presburger/Simplex.h
+++ b/mlir/include/mlir/Analysis/Presburger/Simplex.h
@@ -771,9 +771,10 @@ class Simplex : public SimplexBase {
std::pair<MaybeOptimum<MPInt>, MaybeOptimum<MPInt>>
computeIntegerBounds(ArrayRef<MPInt> coeffs);
- /// Returns false if the given equality has distinct upper and lower bounds
- /// in the simplex.
- bool isValidEquality(const ArrayRef<MPInt> coeffs);
+ /// Check if the simplex takes only one rational value along the
+ /// direction of `coeffs`.
+ /// `this` must be nonempty.
+ bool isFlatAlong(const ArrayRef<MPInt> coeffs);
/// Returns true if the polytope is unbounded, i.e., extends to infinity in
/// some direction. Otherwise, returns false.
diff --git a/mlir/lib/Analysis/Presburger/Barvinok.cpp b/mlir/lib/Analysis/Presburger/Barvinok.cpp
index ae66bedd21cd2..70ec9da4d530f 100644
--- a/mlir/lib/Analysis/Presburger/Barvinok.cpp
+++ b/mlir/lib/Analysis/Presburger/Barvinok.cpp
@@ -161,33 +161,33 @@ std::optional<ParamPoint>
mlir::presburger::detail::solveParametricEquations(FracMatrix equations) {
// equations is a d x (d + p + 1) matrix.
// Each row represents an equation.
- unsigned numEqs = equations.getNumRows();
+ unsigned d = equations.getNumRows();
unsigned numCols = equations.getNumColumns();
// If the determinant is zero, there is no unique solution.
// Thus we return null.
- if (FracMatrix(equations.getSubMatrix(/*fromRow=*/0, /*toRow=*/numEqs - 1,
+ if (FracMatrix(equations.getSubMatrix(/*fromRow=*/0, /*toRow=*/d - 1,
/*fromColumn=*/0,
- /*toColumn=*/numEqs - 1))
+ /*toColumn=*/d - 1))
.determinant() == 0)
return std::nullopt;
- for (unsigned i = 0; i < numEqs; ++i) {
+ for (unsigned i = 0; i < d; ++i) {
// First ensure that the diagonal element is nonzero, by swapping
// it with a row that is non-zero at column i.
- if (equations(i, i) == 0) {
- for (unsigned j = i + 1; j < numEqs; ++j) {
- if (equations(j, i) == 0)
- continue;
- equations.swapRows(j, i);
- break;
- }
+ if (equations(i, i) != 0)
+ continue;
+ for (unsigned j = i + 1; j < d; ++j) {
+ if (equations(j, i) == 0)
+ continue;
+ equations.swapRows(j, i);
+ break;
}
Fraction diagElement = equations(i, i);
// Apply row operations to make all elements except the diagonal to zero.
- for (unsigned j = 0; j < numEqs; ++j) {
+ for (unsigned j = 0; j < d; ++j) {
if (i == j)
continue;
if (equations(j, i) == 0)
@@ -195,12 +195,13 @@ mlir::presburger::detail::solveParametricEquations(FracMatrix equations) {
// Apply row operations to make element (j, i) zero by subtracting the
// ith row, appropriately scaled.
Fraction currentElement = equations(j, i);
- equations.addToRow(j, equations.getRow(i), -currentElement / diagElement);
+ equations.addToRow(/*sourceRow=*/i, /*targetRow=*/j,
+ -currentElement / diagElement);
}
}
// Rescale diagonal elements to 1.
- for (unsigned i = 0; i < numEqs; ++i) {
+ for (unsigned i = 0; i < d; ++i) {
Fraction diagElement = equations(i, i);
for (unsigned j = 0; j < numCols; ++j)
equations(i, j) = equations(i, j) / diagElement;
@@ -215,9 +216,9 @@ mlir::presburger::detail::solveParametricEquations(FracMatrix equations) {
// and so we return the negation of the last p + 1 columns of the matrix.
// We copy these columns and return them.
ParamPoint vertex =
- equations.getSubMatrix(/*fromRow=*/0, /*toRow=*/numEqs - 1,
- /*fromColumn=*/numEqs, /*toColumn=*/numCols - 1);
- for (unsigned i = 0; i < numEqs; ++i)
+ equations.getSubMatrix(/*fromRow=*/0, /*toRow=*/d - 1,
+ /*fromColumn=*/d, /*toColumn=*/numCols - 1);
+ for (unsigned i = 0; i < d; ++i)
vertex.negateRow(i);
return vertex;
@@ -229,7 +230,7 @@ mlir::presburger::detail::solveParametricEquations(FracMatrix equations) {
/// we iterate over the vertex list, each time appending the vertex to the
/// chambers where it is active and creating a new chamber if necessary.
std::vector<std::pair<PresburgerRelation, std::vector<unsigned>>>
-mlir::presburger::detail::chamberDecomposition(
+mlir::presburger::detail::computeChamberDecomposition(
std::vector<PresburgerRelation> activeRegions,
std::vector<ParamPoint> vertices) {
// We maintain a list of regions and their associated vertex sets,
@@ -428,7 +429,7 @@ mlir::presburger::detail::polytopeGeneratingFunction(PolyhedronH poly) {
// they may share "faces" or "edges", but their intersection can only have
// up to numVars-1 dimensions.
std::vector<std::pair<PresburgerRelation, std::vector<unsigned>>> chambers =
- chamberDecomposition(activeRegions, vertices);
+ computeChamberDecomposition(activeRegions, vertices);
// Now, we compute the generating function. For each chamber, we iterate over
// the vertices active in it, and compute the generating function for each
diff --git a/mlir/lib/Analysis/Presburger/IntegerRelation.cpp b/mlir/lib/Analysis/Presburger/IntegerRelation.cpp
index 5895b927435a4..201dc9bd4c8a1 100644
--- a/mlir/lib/Analysis/Presburger/IntegerRelation.cpp
+++ b/mlir/lib/Analysis/Presburger/IntegerRelation.cpp
@@ -26,6 +26,7 @@
#include "llvm/ADT/DenseMap.h"
#include "llvm/ADT/DenseSet.h"
#include "llvm/ADT/STLExtras.h"
+#include "llvm/ADT/Sequence.h"
#include "llvm/ADT/SmallBitVector.h"
#include "llvm/Support/Debug.h"
#include "llvm/Support/raw_ostream.h"
@@ -2499,22 +2500,17 @@ void IntegerRelation::printSpace(raw_ostream &os) const {
}
void IntegerRelation::removeTrivialEqualities() {
- std::vector<bool> isTrivialEquality;
- for (int i = 0, e = getNumEqualities(); i < e; ++i) {
- bool currentIsTrivial =
- llvm::all_of(getEquality(i), [&](MPInt n) -> bool { return (n == 0); });
- isTrivialEquality.push_back(currentIsTrivial);
- }
-
- unsigned pos = 0;
- for (unsigned r = 0, e = getNumEqualities(); r < e; r++)
- if (!isTrivialEquality[r])
- equalities.copyRow(r, pos++);
-
- equalities.resizeVertically(pos);
+ for (int i = getNumEqualities() - 1; i >= 0; --i)
+ if (rangeIsZero(getEquality(i)))
+ removeEquality(i);
}
bool IntegerRelation::isFullDim() {
+ if (getNumVars() == 0)
+ return true;
+ if (isEmpty())
+ return false;
+
// If there is a non-trivial equality, the space cannot be full-dimensional.
removeTrivialEqualities();
if (getNumEqualities() > 0)
@@ -2523,24 +2519,11 @@ bool IntegerRelation::isFullDim() {
// If along the direction of any of the inequalities, the upper and lower
// optima are the same, then the region is not full-dimensional.
Simplex simplex(*this);
- for (unsigned i = 0, e = getNumInequalities(); i < e; i++) {
- auto ineq = inequalities.getRow(i);
- auto upOpt = simplex.computeOptimum(Simplex::Direction::Up, ineq);
- auto downOpt = simplex.computeOptimum(Simplex::Direction::Down, ineq);
-
- if (upOpt.getKind() == OptimumKind::Unbounded ||
- downOpt.getKind() == OptimumKind::Unbounded)
- continue;
-
- // Check if the upper and lower optima are equal.
- if (upOpt.getKind() == OptimumKind::Bounded &&
- downOpt.getKind() == OptimumKind::Bounded && *upOpt == *downOpt)
- return false;
- }
- // If none of the inequalities were such that the upper and lower optima
- // along their direction were equal, then we conclude that the region is full
- // dimensional.
- return true;
+ bool fullDim = llvm::none_of(llvm::seq<int>(0, getNumInequalities()),
+ [&](unsigned i) -> bool {
+ return simplex.isFlatAlong(getInequality(i));
+ });
+ return fullDim;
}
void IntegerRelation::print(raw_ostream &os) const {
diff --git a/mlir/lib/Analysis/Presburger/Matrix.cpp b/mlir/lib/Analysis/Presburger/Matrix.cpp
index 35af156552c34..af4435bda4172 100644
--- a/mlir/lib/Analysis/Presburger/Matrix.cpp
+++ b/mlir/lib/Analysis/Presburger/Matrix.cpp
@@ -358,8 +358,8 @@ template <typename T>
Matrix<T> Matrix<T>::getSubMatrix(unsigned fromRow, unsigned toRow,
unsigned fromColumn,
unsigned toColumn) const {
- assert(toRow >= fromRow && "end of row range must be after beginning!");
- assert(toColumn >= fromColumn && "end of row range must be after beginning!");
+ assert(fromRow <= toRow && "end of row range must be after beginning!");
+ assert(fromColumn <= toColumn && "end of row range must be after beginning!");
Matrix<T> subMatrix(toRow - fromRow + 1, toColumn - fromColumn + 1);
for (unsigned i = fromRow; i <= toRow; i++)
for (unsigned j = fromColumn; j <= toColumn; j++)
diff --git a/mlir/lib/Analysis/Presburger/Simplex.cpp b/mlir/lib/Analysis/Presburger/Simplex.cpp
index f6d24b3f4dbd4..ee3ec87e592c4 100644
--- a/mlir/lib/Analysis/Presburger/Simplex.cpp
+++ b/mlir/lib/Analysis/Presburger/Simplex.cpp
@@ -2104,12 +2104,20 @@ Simplex::computeIntegerBounds(ArrayRef<MPInt> coeffs) {
return {minRoundedUp, maxRoundedDown};
}
-bool Simplex::isValidEquality(ArrayRef<MPInt> coeffs) {
- auto [downOpt, upOpt] = Simplex::computeIntegerBounds(coeffs);
- if (upOpt.getKind() == OptimumKind::Bounded &&
- downOpt.getKind() == OptimumKind::Bounded && *upOpt == *downOpt)
+bool Simplex::isFlatAlong(ArrayRef<MPInt> coeffs) {
+ assert(!isEmpty() && "cannot check for flatness of empty simplex!");
+ auto upOpt = computeOptimum(Simplex::Direction::Up, coeffs);
+ auto downOpt = computeOptimum(Simplex::Direction::Down, coeffs);
+
+ if (upOpt.getKind() != OptimumKind::Bounded)
return false;
- return true;
+ if (downOpt.getKind() != OptimumKind::Bounded)
+ return false;
+
+ // Check if the upper and lower optima are equal.
+ if (*upOpt == *downOpt)
+ return true;
+ return false;
}
void SimplexBase::print(raw_ostream &os) const {
>From 21086534cbb0c9fcc3e691be7983318300076937 Mon Sep 17 00:00:00 2001
From: Abhinav271828 <abhinav.m at research.iiit.ac.in>
Date: Wed, 24 Jan 2024 03:07:08 +0530
Subject: [PATCH 23/41] Minor
---
mlir/lib/Analysis/Presburger/Matrix.cpp | 6 +++---
1 file changed, 3 insertions(+), 3 deletions(-)
diff --git a/mlir/lib/Analysis/Presburger/Matrix.cpp b/mlir/lib/Analysis/Presburger/Matrix.cpp
index af4435bda4172..0ba1dc4107914 100644
--- a/mlir/lib/Analysis/Presburger/Matrix.cpp
+++ b/mlir/lib/Analysis/Presburger/Matrix.cpp
@@ -361,8 +361,8 @@ Matrix<T> Matrix<T>::getSubMatrix(unsigned fromRow, unsigned toRow,
assert(fromRow <= toRow && "end of row range must be after beginning!");
assert(fromColumn <= toColumn && "end of row range must be after beginning!");
Matrix<T> subMatrix(toRow - fromRow + 1, toColumn - fromColumn + 1);
- for (unsigned i = fromRow; i <= toRow; i++)
- for (unsigned j = fromColumn; j <= toColumn; j++)
+ for (unsigned i = fromRow; i <= toRow; ++i)
+ for (unsigned j = fromColumn; j <= toColumn; ++j)
subMatrix(i - fromRow, j - fromColumn) = at(i, j);
return subMatrix;
}
@@ -388,7 +388,7 @@ Matrix<T>::splitByBitset(std::bitset<16> indicator) {
else
rowsForZero.appendExtraRow(getRow(i));
}
- return std::make_pair(rowsForOne, rowsForZero);
+ return {rowsForOne, rowsForZero};
}
template <typename T>
>From 540a0c15b660e7200c535ec15066a1cd25466bb1 Mon Sep 17 00:00:00 2001
From: Abhinav271828 <abhinav.m at research.iiit.ac.in>
Date: Wed, 24 Jan 2024 22:56:20 +0530
Subject: [PATCH 24/41] Minor fixes
---
.../mlir/Analysis/Presburger/Barvinok.h | 5 +-
.../Analysis/Presburger/IntegerRelation.h | 4 +-
.../include/mlir/Analysis/Presburger/Matrix.h | 2 +-
.../Analysis/Presburger/PresburgerRelation.h | 4 ++
.../mlir/Analysis/Presburger/Simplex.h | 3 +-
mlir/lib/Analysis/Presburger/Barvinok.cpp | 51 ++++++++++---------
.../Analysis/Presburger/IntegerRelation.cpp | 8 ++-
mlir/lib/Analysis/Presburger/Matrix.cpp | 5 +-
.../Presburger/PresburgerRelation.cpp | 7 +++
mlir/lib/Analysis/Presburger/Simplex.cpp | 8 ++-
10 files changed, 57 insertions(+), 40 deletions(-)
diff --git a/mlir/include/mlir/Analysis/Presburger/Barvinok.h b/mlir/include/mlir/Analysis/Presburger/Barvinok.h
index 4d3f643187247..f322842e9aff3 100644
--- a/mlir/include/mlir/Analysis/Presburger/Barvinok.h
+++ b/mlir/include/mlir/Analysis/Presburger/Barvinok.h
@@ -101,11 +101,12 @@ std::optional<ParamPoint> solveParametricEquations(FracMatrix equations);
/// In the return type, the vertices are stored by their index in the
/// `vertices` argument, i.e., the set {2, 3, 4} represents the vertex set
/// {vertices[2], vertices[3], vertices[4]}.
+/// The ith relation corresponds to the activity region of the ith vertex set.
/// Note that here, by disjoint, we mean that the intersection is not
/// full-dimensional.
std::vector<std::pair<PresburgerRelation, std::vector<unsigned>>>
-computeChamberDecomposition(std::vector<PresburgerRelation> activeRegions,
- std::vector<ParamPoint> vertices);
+computeChamberDecomposition(ArrayRef<PresburgerRelation> activeRegions,
+ ArrayRef<ParamPoint> vertices);
/// Compute the generating function corresponding to a polytope.
/// All tangent cones of the polytope must be unimodular.
diff --git a/mlir/include/mlir/Analysis/Presburger/IntegerRelation.h b/mlir/include/mlir/Analysis/Presburger/IntegerRelation.h
index 3e76e4bd2cc0c..4c4576bf6d0d3 100644
--- a/mlir/include/mlir/Analysis/Presburger/IntegerRelation.h
+++ b/mlir/include/mlir/Analysis/Presburger/IntegerRelation.h
@@ -716,8 +716,10 @@ class IntegerRelation {
// Verify whether the relation is full-dimensional, i.e.,
// no equality holds for the relation.
+ //
// If there are no variables, it always returns true.
- // If there is at least one variable and the space is empty, it returns false.
+ // If there is at least one variable and the relation is empty, it returns
+ // false.
bool isFullDim();
void print(raw_ostream &os) const;
diff --git a/mlir/include/mlir/Analysis/Presburger/Matrix.h b/mlir/include/mlir/Analysis/Presburger/Matrix.h
index ae781a93a036f..9cd8cebeb2c84 100644
--- a/mlir/include/mlir/Analysis/Presburger/Matrix.h
+++ b/mlir/include/mlir/Analysis/Presburger/Matrix.h
@@ -311,7 +311,7 @@ class FracMatrix : public Matrix<Fraction> {
// calls `y`, usually 3/4.
void LLL(Fraction delta);
- // Multiply each row of the matrix by the LCD of the denominators, thereby
+ // Multiply each row of the matrix by the LCM of the denominators, thereby
// converting it to an integer matrix.
IntMatrix normalizeRows();
};
diff --git a/mlir/include/mlir/Analysis/Presburger/PresburgerRelation.h b/mlir/include/mlir/Analysis/Presburger/PresburgerRelation.h
index c6b00eca90733..35c89bc6e6846 100644
--- a/mlir/include/mlir/Analysis/Presburger/PresburgerRelation.h
+++ b/mlir/include/mlir/Analysis/Presburger/PresburgerRelation.h
@@ -217,6 +217,10 @@ class PresburgerRelation {
/// redundencies.
PresburgerRelation simplify() const;
+ /// Return true if any of the disjuncts is full-dimensional
+ /// (see IntegerRelation::isFullDim()).
+ bool isFullDim() const;
+
/// Print the set's internal state.
void print(raw_ostream &os) const;
void dump() const;
diff --git a/mlir/include/mlir/Analysis/Presburger/Simplex.h b/mlir/include/mlir/Analysis/Presburger/Simplex.h
index 8be9942a3dab9..7ee74c150867c 100644
--- a/mlir/include/mlir/Analysis/Presburger/Simplex.h
+++ b/mlir/include/mlir/Analysis/Presburger/Simplex.h
@@ -773,8 +773,9 @@ class Simplex : public SimplexBase {
/// Check if the simplex takes only one rational value along the
/// direction of `coeffs`.
+ ///
/// `this` must be nonempty.
- bool isFlatAlong(const ArrayRef<MPInt> coeffs);
+ bool isFlatAlong(ArrayRef<MPInt> coeffs);
/// Returns true if the polytope is unbounded, i.e., extends to infinity in
/// some direction. Otherwise, returns false.
diff --git a/mlir/lib/Analysis/Presburger/Barvinok.cpp b/mlir/lib/Analysis/Presburger/Barvinok.cpp
index 70ec9da4d530f..1d22880393846 100644
--- a/mlir/lib/Analysis/Presburger/Barvinok.cpp
+++ b/mlir/lib/Analysis/Presburger/Barvinok.cpp
@@ -172,6 +172,8 @@ mlir::presburger::detail::solveParametricEquations(FracMatrix equations) {
.determinant() == 0)
return std::nullopt;
+ // Perform row operations to make each column all zeros except for the
+ // diagonal element, which is made to be one.
for (unsigned i = 0; i < d; ++i) {
// First ensure that the diagonal element is nonzero, by swapping
// it with a row that is non-zero at column i.
@@ -196,7 +198,7 @@ mlir::presburger::detail::solveParametricEquations(FracMatrix equations) {
// ith row, appropriately scaled.
Fraction currentElement = equations(j, i);
equations.addToRow(/*sourceRow=*/i, /*targetRow=*/j,
- -currentElement / diagElement);
+ /*scale=*/-currentElement / diagElement);
}
}
@@ -229,10 +231,14 @@ mlir::presburger::detail::solveParametricEquations(FracMatrix equations) {
/// We maintain a list of pairwise disjoint chambers and their vertex-sets;
/// we iterate over the vertex list, each time appending the vertex to the
/// chambers where it is active and creating a new chamber if necessary.
+///
+/// Given the region each vertex is active in, for each subset of vertices,
+/// the region that precisely this subset is in, is the intersection of the
+/// regions that these are active in, intersected with the complements of the
+/// remaining regions.
std::vector<std::pair<PresburgerRelation, std::vector<unsigned>>>
mlir::presburger::detail::computeChamberDecomposition(
- std::vector<PresburgerRelation> activeRegions,
- std::vector<ParamPoint> vertices) {
+ ArrayRef<PresburgerRelation> activeRegions, ArrayRef<ParamPoint> vertices) {
// We maintain a list of regions and their associated vertex sets,
// initialized with the first vertex and its corresponding activity region.
std::vector<std::pair<PresburgerRelation, std::vector<unsigned>>> chambers = {
@@ -241,8 +247,6 @@ mlir::presburger::detail::computeChamberDecomposition(
// of indices. Thus the set {2, 3, 4} represents the vertex set
// {vertices[2], vertices[3], vertices[4]}.
- std::vector<std::pair<PresburgerRelation, std::vector<unsigned>>> newChambers;
-
// We iterate over the vertex set.
// For each vertex v_j and its activity region R_j,
// we examine all the current chambers R_i.
@@ -253,37 +257,35 @@ mlir::presburger::detail::computeChamberDecomposition(
// Once we have examined all R_i, we add a final chamber
// R_j - (union of all existing chambers),
// in which only v_j is active.
+
+ // At each step, we define a new chamber list after considering vertex v_j,
+ // replacing and appending chambers as discussed above.
+ std::vector<std::pair<PresburgerRelation, std::vector<unsigned>>> newChambers;
for (unsigned j = 1, e = vertices.size(); j < e; j++) {
newChambers.clear();
PresburgerRelation newRegion = activeRegions[j];
- ParamPoint newVertex = vertices[j];
-
- for (unsigned i = 0, f = chambers.size(); i < f; i++) {
- auto [currentRegion, currentVertices] = chambers[i];
+ for (auto [currentRegion, currentVertices] : chambers) {
// First, we check if the intersection of R_j and R_i.
// It is a disjoint union of convex regions in the parameter space,
// and so we know that it is full-dimensional if any of the disjuncts
// is full-dimensional.
PresburgerRelation intersection = currentRegion.intersect(newRegion);
- bool isFullDim = intersection.getNumRangeVars() == 0 ||
- llvm::any_of(intersection.getAllDisjuncts(),
- [&](IntegerRelation disjunct) -> bool {
- return disjunct.isFullDim();
- });
+ bool isFullDim =
+ intersection.getNumRangeVars() == 0 || intersection.isFullDim();
// If the intersection is not full-dimensional, we do not modify
// the chamber list.
if (!isFullDim)
- newChambers.push_back(chambers[i]);
+ newChambers.emplace_back(currentRegion, currentVertices);
else {
// If it is, we add the intersection and the difference as new chambers.
PresburgerRelation subtraction = currentRegion.subtract(newRegion);
- newChambers.push_back(std::make_pair(subtraction, currentVertices));
+ newChambers.emplace_back(subtraction, currentVertices);
currentVertices.push_back(j);
- newChambers.push_back(std::make_pair(intersection, currentVertices));
+ newChambers.emplace_back(intersection, currentVertices);
}
}
@@ -441,27 +443,28 @@ mlir::presburger::detail::polytopeGeneratingFunction(PolyhedronH poly) {
auto [currentRegion, currentVertices] = chamber;
GeneratingFunction chamberGf(numSymbols, {}, {}, {});
for (unsigned i : currentVertices) {
- // We collect the inequalities corresponding to each vertex.
+ // We collect the inequalities corresponding to each vertex to compute
+ // the tangent cone at that vertex.
// We only need the coefficients of the variables (NOT the parameters)
// as the generating function only depends on these.
// We translate the cones to be pointed at the origin by making the
// constant terms zero.
- ConeH tgtCone = defineHRep(numVars);
- for (unsigned j = 0; j < numVars; j++) {
- for (unsigned k = 0; k < numVars; k++)
+ ConeH tangentCone = defineHRep(numVars);
+ for (unsigned j = 0; j < numVars; ++j) {
+ for (unsigned k = 0; k < numVars; ++k)
ineq[k] = subsets[i](j, k);
- tgtCone.addInequality(ineq);
+ tangentCone.addInequality(ineq);
}
// We assume that the tangent cone is unimodular.
SmallVector<std::pair<int, ConeH>, 4> unimodCones = {
- std::make_pair(1, tgtCone)};
+ std::make_pair(1, tangentCone)};
for (std::pair<int, ConeH> signedCone : unimodCones) {
auto [sign, cone] = signedCone;
chamberGf = chamberGf +
unimodularConeGeneratingFunction(vertices[i], sign, cone);
}
}
- gf.push_back(std::make_pair(currentRegion, chamberGf));
+ gf.emplace_back(currentRegion, chamberGf);
}
return gf;
}
diff --git a/mlir/lib/Analysis/Presburger/IntegerRelation.cpp b/mlir/lib/Analysis/Presburger/IntegerRelation.cpp
index 201dc9bd4c8a1..ab50f0dbec2cf 100644
--- a/mlir/lib/Analysis/Presburger/IntegerRelation.cpp
+++ b/mlir/lib/Analysis/Presburger/IntegerRelation.cpp
@@ -2519,11 +2519,9 @@ bool IntegerRelation::isFullDim() {
// If along the direction of any of the inequalities, the upper and lower
// optima are the same, then the region is not full-dimensional.
Simplex simplex(*this);
- bool fullDim = llvm::none_of(llvm::seq<int>(0, getNumInequalities()),
- [&](unsigned i) -> bool {
- return simplex.isFlatAlong(getInequality(i));
- });
- return fullDim;
+ return llvm::none_of(llvm::seq<int>(0, getNumInequalities()), [&](int i) {
+ return simplex.isFlatAlong(getInequality(i));
+ });
}
void IntegerRelation::print(raw_ostream &os) const {
diff --git a/mlir/lib/Analysis/Presburger/Matrix.cpp b/mlir/lib/Analysis/Presburger/Matrix.cpp
index 0ba1dc4107914..6780288463a7d 100644
--- a/mlir/lib/Analysis/Presburger/Matrix.cpp
+++ b/mlir/lib/Analysis/Presburger/Matrix.cpp
@@ -359,7 +359,10 @@ Matrix<T> Matrix<T>::getSubMatrix(unsigned fromRow, unsigned toRow,
unsigned fromColumn,
unsigned toColumn) const {
assert(fromRow <= toRow && "end of row range must be after beginning!");
- assert(fromColumn <= toColumn && "end of row range must be after beginning!");
+ assert(toRow < nRows && "end of row range out of bounds!");
+ assert(fromColumn <= toColumn &&
+ "end of column range must be after beginning!");
+ assert(toColumn < nColumns && "end of column range out of bounds!");
Matrix<T> subMatrix(toRow - fromRow + 1, toColumn - fromColumn + 1);
for (unsigned i = fromRow; i <= toRow; ++i)
for (unsigned j = fromColumn; j <= toColumn; ++j)
diff --git a/mlir/lib/Analysis/Presburger/PresburgerRelation.cpp b/mlir/lib/Analysis/Presburger/PresburgerRelation.cpp
index 787fc1c659a12..ee56ce696e86e 100644
--- a/mlir/lib/Analysis/Presburger/PresburgerRelation.cpp
+++ b/mlir/lib/Analysis/Presburger/PresburgerRelation.cpp
@@ -1041,6 +1041,13 @@ PresburgerRelation PresburgerRelation::simplify() const {
return result;
}
+bool PresburgerRelation::isFullDim() const {
+ bool fullDim = llvm::any_of(getAllDisjuncts(), [&](IntegerRelation disjunct) {
+ return disjunct.isFullDim();
+ });
+ return fullDim;
+}
+
void PresburgerRelation::print(raw_ostream &os) const {
os << "Number of Disjuncts: " << getNumDisjuncts() << "\n";
for (const IntegerRelation &disjunct : disjuncts) {
diff --git a/mlir/lib/Analysis/Presburger/Simplex.cpp b/mlir/lib/Analysis/Presburger/Simplex.cpp
index ee3ec87e592c4..592d96ee766fb 100644
--- a/mlir/lib/Analysis/Presburger/Simplex.cpp
+++ b/mlir/lib/Analysis/Presburger/Simplex.cpp
@@ -2109,15 +2109,13 @@ bool Simplex::isFlatAlong(ArrayRef<MPInt> coeffs) {
auto upOpt = computeOptimum(Simplex::Direction::Up, coeffs);
auto downOpt = computeOptimum(Simplex::Direction::Down, coeffs);
- if (upOpt.getKind() != OptimumKind::Bounded)
+ if (!upOpt.isBounded())
return false;
- if (downOpt.getKind() != OptimumKind::Bounded)
+ if (!downOpt.isBounded())
return false;
// Check if the upper and lower optima are equal.
- if (*upOpt == *downOpt)
- return true;
- return false;
+ return *upOpt == *downOpt;
}
void SimplexBase::print(raw_ostream &os) const {
>From 5dc100f478a14635a41f11ac827ed7212651ca38 Mon Sep 17 00:00:00 2001
From: Abhinav271828 <abhinav.m at research.iiit.ac.in>
Date: Wed, 24 Jan 2024 22:58:26 +0530
Subject: [PATCH 25/41] Use const
---
mlir/include/mlir/Analysis/Presburger/Matrix.h | 2 +-
mlir/lib/Analysis/Presburger/Matrix.cpp | 2 +-
2 files changed, 2 insertions(+), 2 deletions(-)
diff --git a/mlir/include/mlir/Analysis/Presburger/Matrix.h b/mlir/include/mlir/Analysis/Presburger/Matrix.h
index 9cd8cebeb2c84..fe260568a22d8 100644
--- a/mlir/include/mlir/Analysis/Presburger/Matrix.h
+++ b/mlir/include/mlir/Analysis/Presburger/Matrix.h
@@ -313,7 +313,7 @@ class FracMatrix : public Matrix<Fraction> {
// Multiply each row of the matrix by the LCM of the denominators, thereby
// converting it to an integer matrix.
- IntMatrix normalizeRows();
+ IntMatrix normalizeRows() const;
};
} // namespace presburger
diff --git a/mlir/lib/Analysis/Presburger/Matrix.cpp b/mlir/lib/Analysis/Presburger/Matrix.cpp
index 6780288463a7d..6d734c232f182 100644
--- a/mlir/lib/Analysis/Presburger/Matrix.cpp
+++ b/mlir/lib/Analysis/Presburger/Matrix.cpp
@@ -729,7 +729,7 @@ void FracMatrix::LLL(Fraction delta) {
}
}
-IntMatrix FracMatrix::normalizeRows() {
+IntMatrix FracMatrix::normalizeRows() const {
unsigned numRows = getNumRows();
unsigned numColumns = getNumColumns();
IntMatrix normalized(numRows, numColumns);
>From aaedcc0c149cce62917f5d9eaa1626041f57fea3 Mon Sep 17 00:00:00 2001
From: Abhinav271828 <abhinav.m at research.iiit.ac.in>
Date: Wed, 24 Jan 2024 23:01:12 +0530
Subject: [PATCH 26/41] Rename
---
mlir/include/mlir/Analysis/Presburger/Barvinok.h | 7 ++++---
mlir/lib/Analysis/Presburger/Barvinok.cpp | 9 +++++----
mlir/unittests/Analysis/Presburger/BarvinokTest.cpp | 9 +++++----
3 files changed, 14 insertions(+), 11 deletions(-)
diff --git a/mlir/include/mlir/Analysis/Presburger/Barvinok.h b/mlir/include/mlir/Analysis/Presburger/Barvinok.h
index f322842e9aff3..729a49b0e8e32 100644
--- a/mlir/include/mlir/Analysis/Presburger/Barvinok.h
+++ b/mlir/include/mlir/Analysis/Presburger/Barvinok.h
@@ -85,8 +85,9 @@ ConeH getDual(ConeV cone);
/// Compute the generating function for a unimodular cone.
/// The input cone must be unimodular; it assert-fails otherwise.
-GeneratingFunction unimodularConeGeneratingFunction(ParamPoint vertex, int sign,
- ConeH cone);
+GeneratingFunction computeUnimodularConeGeneratingFunction(ParamPoint vertex,
+ int sign,
+ ConeH cone);
/// Find the solution of a set of equations that express affine constraints
/// between a set of variables and a set of parameters. The solution expresses
@@ -111,7 +112,7 @@ computeChamberDecomposition(ArrayRef<PresburgerRelation> activeRegions,
/// Compute the generating function corresponding to a polytope.
/// All tangent cones of the polytope must be unimodular.
std::vector<std::pair<PresburgerRelation, GeneratingFunction>>
-polytopeGeneratingFunction(PolyhedronH poly);
+computePolytopeGeneratingFunction(PolyhedronH poly);
/// Find a vector that is not orthogonal to any of the given vectors,
/// i.e., has nonzero dot product with those of the given vectors
diff --git a/mlir/lib/Analysis/Presburger/Barvinok.cpp b/mlir/lib/Analysis/Presburger/Barvinok.cpp
index 1d22880393846..b997b48d960e3 100644
--- a/mlir/lib/Analysis/Presburger/Barvinok.cpp
+++ b/mlir/lib/Analysis/Presburger/Barvinok.cpp
@@ -77,7 +77,8 @@ MPInt mlir::presburger::detail::getIndex(ConeV cone) {
/// num is computed by expressing the vertex as a weighted
/// sum of the generators, and then taking the floor of the
/// coefficients.
-GeneratingFunction mlir::presburger::detail::unimodularConeGeneratingFunction(
+GeneratingFunction
+mlir::presburger::detail::computeUnimodularConeGeneratingFunction(
ParamPoint vertex, int sign, ConeH cone) {
// Consider a cone with H-representation [0 -1].
// [-1 -2]
@@ -328,7 +329,7 @@ mlir::presburger::detail::computeChamberDecomposition(
/// polytopes using Barvinok's rational functions." Algorithmica 48 (2007):
/// 37-66.
std::vector<std::pair<PresburgerRelation, GeneratingFunction>>
-mlir::presburger::detail::polytopeGeneratingFunction(PolyhedronH poly) {
+mlir::presburger::detail::computePolytopeGeneratingFunction(PolyhedronH poly) {
unsigned numVars = poly.getNumRangeVars();
unsigned numSymbols = poly.getNumSymbolVars();
unsigned numIneqs = poly.getNumInequalities();
@@ -460,8 +461,8 @@ mlir::presburger::detail::polytopeGeneratingFunction(PolyhedronH poly) {
std::make_pair(1, tangentCone)};
for (std::pair<int, ConeH> signedCone : unimodCones) {
auto [sign, cone] = signedCone;
- chamberGf = chamberGf +
- unimodularConeGeneratingFunction(vertices[i], sign, cone);
+ chamberGf = chamberGf + computeUnimodularConeGeneratingFunction(
+ vertices[i], sign, cone);
}
}
gf.emplace_back(currentRegion, chamberGf);
diff --git a/mlir/unittests/Analysis/Presburger/BarvinokTest.cpp b/mlir/unittests/Analysis/Presburger/BarvinokTest.cpp
index 09301d5e54e0e..50a6cfa0fb6cc 100644
--- a/mlir/unittests/Analysis/Presburger/BarvinokTest.cpp
+++ b/mlir/unittests/Analysis/Presburger/BarvinokTest.cpp
@@ -59,7 +59,8 @@ TEST(BarvinokTest, unimodularConeGeneratingFunction) {
ParamPoint vertex =
makeFracMatrix(2, 3, {{2, 2, 0}, {-1, -Fraction(1, 2), 1}});
- GeneratingFunction gf = unimodularConeGeneratingFunction(vertex, 1, cone);
+ GeneratingFunction gf =
+ computeUnimodularConeGeneratingFunction(vertex, 1, cone);
EXPECT_EQ_REPR_GENERATINGFUNCTION(
gf, GeneratingFunction(
@@ -74,7 +75,7 @@ TEST(BarvinokTest, unimodularConeGeneratingFunction) {
vertex = makeFracMatrix(3, 2, {{5, 2}, {6, 2}, {7, 1}});
- gf = unimodularConeGeneratingFunction(vertex, 1, cone);
+ gf = computeUnimodularConeGeneratingFunction(vertex, 1, cone);
EXPECT_EQ_REPR_GENERATINGFUNCTION(
gf,
@@ -252,7 +253,7 @@ TEST(BarvinokTest, computeNumTermsPolytope) {
poly.addInequality(ineqs.getRow(i));
std::vector<std::pair<PresburgerRelation, GeneratingFunction>> count =
- polytopeGeneratingFunction(poly);
+ computePolytopeGeneratingFunction(poly);
// There is only one chamber, as it is non-parametric.
EXPECT_EQ(count.size(), 1u);
@@ -281,7 +282,7 @@ TEST(BarvinokTest, computeNumTermsPolytope) {
for (unsigned i = 0; i < 3; i++)
poly.addInequality(ineqs.getRow(i));
- count = polytopeGeneratingFunction(poly);
+ count = computePolytopeGeneratingFunction(poly);
// There is only one chamber: p ≥ 0
EXPECT_EQ(count.size(), 1u);
>From 4f9e7f3c5b7c6435a87af787e91715255fc19b75 Mon Sep 17 00:00:00 2001
From: Abhinav271828 <abhinav.m at research.iiit.ac.in>
Date: Wed, 24 Jan 2024 23:42:53 +0530
Subject: [PATCH 27/41] shift decl
---
mlir/lib/Analysis/Presburger/Barvinok.cpp | 6 +-----
1 file changed, 1 insertion(+), 5 deletions(-)
diff --git a/mlir/lib/Analysis/Presburger/Barvinok.cpp b/mlir/lib/Analysis/Presburger/Barvinok.cpp
index b997b48d960e3..7459692763ba3 100644
--- a/mlir/lib/Analysis/Presburger/Barvinok.cpp
+++ b/mlir/lib/Analysis/Presburger/Barvinok.cpp
@@ -338,11 +338,6 @@ mlir::presburger::detail::computePolytopeGeneratingFunction(PolyhedronH poly) {
// 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 x p +
- // activeRegionConstant ≥ 0. The active region is a polyhedron in parameter
- // space.
- FracMatrix activeRegion(numIneqs - numVars, numSymbols + 1);
-
// These vectors store lists of
// subsets of inequalities,
// the vertices corresponding to them, and
@@ -411,6 +406,7 @@ mlir::presburger::detail::computePolytopeGeneratingFunction(PolyhedronH poly) {
// 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].
+ FracMatrix activeRegion(numIneqs - numVars, numSymbols + 1);
for (unsigned i = 0; i < numIneqs - numVars; i++) {
activeRegion.setRow(i, (*vertex).preMultiplyWithRow(a2.getRow(i)));
activeRegion.addToRow(i, b2c2.getRow(i), 1);
>From cceb81ac976615d2273df187d20ed8ab03feabb4 Mon Sep 17 00:00:00 2001
From: Abhinav271828 <abhinav.m at research.iiit.ac.in>
Date: Thu, 25 Jan 2024 01:58:21 +0530
Subject: [PATCH 28/41] Use PresburgerSet, refactor chamberDecomp, add test
---
.../mlir/Analysis/Presburger/Barvinok.h | 8 +-
mlir/lib/Analysis/Presburger/Barvinok.cpp | 90 ++++++++-----------
.../Analysis/Presburger/BarvinokTest.cpp | 45 +++++++++-
3 files changed, 83 insertions(+), 60 deletions(-)
diff --git a/mlir/include/mlir/Analysis/Presburger/Barvinok.h b/mlir/include/mlir/Analysis/Presburger/Barvinok.h
index 729a49b0e8e32..389be6914f3cc 100644
--- a/mlir/include/mlir/Analysis/Presburger/Barvinok.h
+++ b/mlir/include/mlir/Analysis/Presburger/Barvinok.h
@@ -95,7 +95,7 @@ GeneratingFunction computeUnimodularConeGeneratingFunction(ParamPoint vertex,
/// If there is no solution, return null.
std::optional<ParamPoint> solveParametricEquations(FracMatrix equations);
-/// Given a list of possibly intersecting regions (PresburgerRelations) and
+/// Given a list of possibly intersecting regions (PresburgerSet) 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.
@@ -105,13 +105,13 @@ std::optional<ParamPoint> solveParametricEquations(FracMatrix equations);
/// The ith relation corresponds to the activity region of the ith vertex set.
/// Note that here, by disjoint, we mean that the intersection is not
/// full-dimensional.
-std::vector<std::pair<PresburgerRelation, std::vector<unsigned>>>
-computeChamberDecomposition(ArrayRef<PresburgerRelation> activeRegions,
+std::vector<std::pair<PresburgerSet, std::vector<unsigned>>>
+computeChamberDecomposition(ArrayRef<PresburgerSet> activeRegions,
ArrayRef<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>>
+std::vector<std::pair<PresburgerSet, GeneratingFunction>>
computePolytopeGeneratingFunction(PolyhedronH poly);
/// Find a vector that is not orthogonal to any of the given vectors,
diff --git a/mlir/lib/Analysis/Presburger/Barvinok.cpp b/mlir/lib/Analysis/Presburger/Barvinok.cpp
index 7459692763ba3..329a983ecf8c6 100644
--- a/mlir/lib/Analysis/Presburger/Barvinok.cpp
+++ b/mlir/lib/Analysis/Presburger/Barvinok.cpp
@@ -231,19 +231,20 @@ mlir::presburger::detail::solveParametricEquations(FracMatrix equations) {
/// 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.
+/// chambers where it is active and separating the chambers where it is not.
///
/// Given the region each vertex is active in, for each subset of vertices,
/// the region that precisely this subset is in, is the intersection of the
/// regions that these are active in, intersected with the complements of the
/// remaining regions.
-std::vector<std::pair<PresburgerRelation, std::vector<unsigned>>>
+std::vector<std::pair<PresburgerSet, std::vector<unsigned>>>
mlir::presburger::detail::computeChamberDecomposition(
- ArrayRef<PresburgerRelation> activeRegions, ArrayRef<ParamPoint> vertices) {
+ ArrayRef<PresburgerSet> activeRegions, ArrayRef<ParamPoint> vertices) {
// We maintain a list of regions and their associated vertex sets,
- // initialized with the first vertex and its corresponding activity region.
- std::vector<std::pair<PresburgerRelation, std::vector<unsigned>>> chambers = {
- std::make_pair(activeRegions[0], std::vector({0u}))};
+ // initialized with the universe and the empty list of vertices.
+ std::vector<std::pair<PresburgerSet, std::vector<unsigned>>> chambers = {
+ {PresburgerSet::getUniverse(activeRegions[0].getSpace()),
+ std::vector<unsigned>({})}};
// 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]}.
@@ -255,59 +256,39 @@ mlir::presburger::detail::computeChamberDecomposition(
// 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.
// At each step, we define a new chamber list after considering vertex v_j,
// replacing and appending chambers as discussed above.
- std::vector<std::pair<PresburgerRelation, std::vector<unsigned>>> newChambers;
- for (unsigned j = 1, e = vertices.size(); j < e; j++) {
- newChambers.clear();
-
- PresburgerRelation newRegion = activeRegions[j];
+ // The loop has the invariant that the union over all the chambers gives the
+ // universe at every step (modulo lower-dimensional spaces).
+ for (unsigned j = 0, e = vertices.size(); j < e; j++) {
+ std::vector<std::pair<PresburgerSet, std::vector<unsigned>>> newChambers;
for (auto [currentRegion, currentVertices] : chambers) {
// First, we check if the intersection of R_j and R_i.
// It is a disjoint union of convex regions in the parameter space,
// and so we know that it is full-dimensional if any of the disjuncts
// is full-dimensional.
- PresburgerRelation intersection = currentRegion.intersect(newRegion);
- bool isFullDim =
- intersection.getNumRangeVars() == 0 || intersection.isFullDim();
+ PresburgerSet intersection = currentRegion.intersect(activeRegions[j]);
// If the intersection is not full-dimensional, we do not modify
// the chamber list.
- if (!isFullDim)
+ if (!intersection.isFullDim()) {
newChambers.emplace_back(currentRegion, currentVertices);
- else {
- // If it is, we add the intersection and the difference as new chambers.
- PresburgerRelation subtraction = currentRegion.subtract(newRegion);
- newChambers.emplace_back(subtraction, currentVertices);
-
- currentVertices.push_back(j);
- newChambers.emplace_back(intersection, currentVertices);
+ continue;
}
- }
- // 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);
+ // If it is, we add the intersection as a chamber where this vertex is
+ // active. We also add the difference if it is full-dimensional.
+ PresburgerSet difference = currentRegion.subtract(activeRegions[j]);
+ if (difference.isFullDim())
+ newChambers.emplace_back(difference, currentVertices);
+
+ currentVertices.push_back(j);
+ newChambers.emplace_back(intersection, currentVertices);
}
+
+ chambers = std::move(newChambers);
}
return chambers;
@@ -328,7 +309,7 @@ mlir::presburger::detail::computeChamberDecomposition(
/// 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>>
+std::vector<std::pair<PresburgerSet, GeneratingFunction>>
mlir::presburger::detail::computePolytopeGeneratingFunction(PolyhedronH poly) {
unsigned numVars = poly.getNumRangeVars();
unsigned numSymbols = poly.getNumSymbolVars();
@@ -336,7 +317,7 @@ mlir::presburger::detail::computePolytopeGeneratingFunction(PolyhedronH poly) {
// 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({});
+ std::vector<std::pair<PresburgerSet, GeneratingFunction>> gf({});
// These vectors store lists of
// subsets of inequalities,
@@ -344,7 +325,7 @@ mlir::presburger::detail::computePolytopeGeneratingFunction(PolyhedronH poly) {
// the active regions of the vertices, in order.
std::vector<IntMatrix> subsets;
std::vector<ParamPoint> vertices;
- std::vector<PresburgerRelation> activeRegions;
+ std::vector<PresburgerSet> activeRegions;
FracMatrix a2(numIneqs - numVars, numVars);
FracMatrix b2c2(numIneqs - numVars, numSymbols + 1);
@@ -413,13 +394,13 @@ mlir::presburger::detail::computePolytopeGeneratingFunction(PolyhedronH poly) {
}
// We convert the representation of the active region to an integers-only
- // form so as to store it as an PresburgerRelation.
+ // form so as to store it as an PresburgerSet.
IntMatrix activeRegionNorm = activeRegion.normalizeRows();
- IntegerRelation activeRegionRel =
- IntegerRelation(PresburgerSpace::getRelationSpace(0, numSymbols, 0, 0));
+ IntegerPolyhedron activeRegionRel = IntegerPolyhedron(
+ PresburgerSpace::getRelationSpace(0, numSymbols, 0, 0));
for (unsigned i = 0, e = activeRegion.getNumRows(); i < e; ++i)
activeRegionRel.addInequality(activeRegionNorm.getRow(i));
- activeRegions.push_back(PresburgerRelation(activeRegionRel));
+ activeRegions.push_back(PresburgerSet(activeRegionRel));
}
// Now, we use Clauss-Loechner decomposition to identify regions in parameter
@@ -427,7 +408,7 @@ mlir::presburger::detail::computePolytopeGeneratingFunction(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.
- std::vector<std::pair<PresburgerRelation, std::vector<unsigned>>> chambers =
+ std::vector<std::pair<PresburgerSet, std::vector<unsigned>>> chambers =
computeChamberDecomposition(activeRegions, vertices);
// Now, we compute the generating function. For each chamber, we iterate over
@@ -435,7 +416,7 @@ mlir::presburger::detail::computePolytopeGeneratingFunction(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 :
+ for (const std::pair<PresburgerSet, std::vector<unsigned>> &chamber :
chambers) {
auto [currentRegion, currentVertices] = chamber;
GeneratingFunction chamberGf(numSymbols, {}, {}, {});
@@ -452,7 +433,10 @@ mlir::presburger::detail::computePolytopeGeneratingFunction(PolyhedronH poly) {
ineq[k] = subsets[i](j, k);
tangentCone.addInequality(ineq);
}
- // We assume that the tangent cone is unimodular.
+ // We assume that the tangent cone is unimodular, so there is no need
+ // to decompose it.
+ // In the general case, the unimodular decomposition may have several
+ // cones.
SmallVector<std::pair<int, ConeH>, 4> unimodCones = {
std::make_pair(1, tangentCone)};
for (std::pair<int, ConeH> signedCone : unimodCones) {
diff --git a/mlir/unittests/Analysis/Presburger/BarvinokTest.cpp b/mlir/unittests/Analysis/Presburger/BarvinokTest.cpp
index 50a6cfa0fb6cc..2ab39d0e10ea5 100644
--- a/mlir/unittests/Analysis/Presburger/BarvinokTest.cpp
+++ b/mlir/unittests/Analysis/Presburger/BarvinokTest.cpp
@@ -252,7 +252,7 @@ TEST(BarvinokTest, computeNumTermsPolytope) {
for (unsigned i = 0; i < 6; i++)
poly.addInequality(ineqs.getRow(i));
- std::vector<std::pair<PresburgerRelation, GeneratingFunction>> count =
+ std::vector<std::pair<PresburgerSet, GeneratingFunction>> count =
computePolytopeGeneratingFunction(poly);
// There is only one chamber, as it is non-parametric.
EXPECT_EQ(count.size(), 1u);
@@ -284,9 +284,9 @@ TEST(BarvinokTest, computeNumTermsPolytope) {
count = computePolytopeGeneratingFunction(poly);
// There is only one chamber: p ≥ 0
- EXPECT_EQ(count.size(), 1u);
+ EXPECT_EQ(count.size(), 2u);
- gf = count[0].second;
+ gf = count[1].second;
EXPECT_EQ_REPR_GENERATINGFUNCTION(
gf, GeneratingFunction(
1, {1, 1, 1},
@@ -294,4 +294,43 @@ TEST(BarvinokTest, computeNumTermsPolytope) {
makeFracMatrix(2, 2, {{0, 1}, {0, 0}}),
makeFracMatrix(2, 2, {{0, 0}, {0, 0}})},
{{{-1, 1}, {-1, 0}}, {{1, -1}, {0, -1}}, {{1, 0}, {0, 1}}}));
+
+ ineqs = makeIntMatrix(6, 5,
+ {{1, 0, 0, 1, 0},
+ {0, 1, 0, 1, 0},
+ {0, 0, 1, 1, 0},
+ {-1, 0, 0, 1, 0},
+ {0, -1, 0, 1, 0},
+ {0, 0, -1, 1, 0}});
+ poly = defineHRep(3, 1);
+ for (unsigned i = 0; i < 6; i++)
+ poly.addInequality(ineqs.getRow(i));
+
+ count = computePolytopeGeneratingFunction(poly);
+
+ EXPECT_EQ(count.size(), 2u);
+
+ gf = count[0].second;
+
+ // Cartesian product of a cube with side M and a right triangle with side N.
+ ineqs = makeIntMatrix(9, 8,
+ {{1, 0, 0, 0, 0, 1, 0, 0},
+ {0, 1, 0, 0, 0, 1, 0, 0},
+ {0, 0, 1, 0, 0, 1, 0, 0},
+ {-1, 0, 0, 0, 0, 1, 0, 0},
+ {0, -1, 0, 0, 0, 1, 0, 0},
+ {0, 0, -1, 0, 0, 1, 0, 0},
+ {0, 0, 0, 1, 0, 0, 0, 0},
+ {0, 0, 0, 0, 1, 0, 0, 0},
+ {0, 0, 0, -1, -1, 0, 1, 0}});
+ poly = defineHRep(5, 2);
+ for (unsigned i = 0; i < 9; i++)
+ poly.addInequality(ineqs.getRow(i));
+
+ count = computePolytopeGeneratingFunction(poly);
+
+ EXPECT_EQ(count.size(), 2u);
+
+ gf = count[1].second;
+ EXPECT_EQ(gf.getNumerators().size(), 24u);
}
>From ac60d0f217fa0817f0e2634a450c99f75eb6eea2 Mon Sep 17 00:00:00 2001
From: Abhinav271828 <abhinav.m at research.iiit.ac.in>
Date: Thu, 25 Jan 2024 12:34:29 +0530
Subject: [PATCH 29/41] Minor fixes
---
.../mlir/Analysis/Presburger/PresburgerRelation.h | 5 +++--
mlir/lib/Analysis/Presburger/Barvinok.cpp | 9 +++------
mlir/lib/Analysis/Presburger/IntegerRelation.cpp | 4 ++--
mlir/lib/Analysis/Presburger/PresburgerRelation.cpp | 3 +--
mlir/lib/Analysis/Presburger/Simplex.cpp | 1 -
5 files changed, 9 insertions(+), 13 deletions(-)
diff --git a/mlir/include/mlir/Analysis/Presburger/PresburgerRelation.h b/mlir/include/mlir/Analysis/Presburger/PresburgerRelation.h
index 35c89bc6e6846..7d28176f477a3 100644
--- a/mlir/include/mlir/Analysis/Presburger/PresburgerRelation.h
+++ b/mlir/include/mlir/Analysis/Presburger/PresburgerRelation.h
@@ -217,8 +217,9 @@ class PresburgerRelation {
/// redundencies.
PresburgerRelation simplify() const;
- /// Return true if any of the disjuncts is full-dimensional
- /// (see IntegerRelation::isFullDim()).
+ /// Return true if any of the disjuncts is full-dimensional.
+ /// We find if a disjunct is full-dimensional by checking if it is flat
+ /// along the dimension of any of its inequalities.
bool isFullDim() const;
/// Print the set's internal state.
diff --git a/mlir/lib/Analysis/Presburger/Barvinok.cpp b/mlir/lib/Analysis/Presburger/Barvinok.cpp
index 329a983ecf8c6..e289f3884f593 100644
--- a/mlir/lib/Analysis/Presburger/Barvinok.cpp
+++ b/mlir/lib/Analysis/Presburger/Barvinok.cpp
@@ -266,9 +266,6 @@ mlir::presburger::detail::computeChamberDecomposition(
for (auto [currentRegion, currentVertices] : chambers) {
// First, we check if the intersection of R_j and R_i.
- // It is a disjoint union of convex regions in the parameter space,
- // and so we know that it is full-dimensional if any of the disjuncts
- // is full-dimensional.
PresburgerSet intersection = currentRegion.intersect(activeRegions[j]);
// If the intersection is not full-dimensional, we do not modify
@@ -389,12 +386,12 @@ mlir::presburger::detail::computePolytopeGeneratingFunction(PolyhedronH poly) {
// and add each row of [B2 | c2].
FracMatrix activeRegion(numIneqs - numVars, numSymbols + 1);
for (unsigned i = 0; i < numIneqs - numVars; i++) {
- activeRegion.setRow(i, (*vertex).preMultiplyWithRow(a2.getRow(i)));
+ activeRegion.setRow(i, vertex->preMultiplyWithRow(a2.getRow(i)));
activeRegion.addToRow(i, b2c2.getRow(i), 1);
}
// We convert the representation of the active region to an integers-only
- // form so as to store it as an PresburgerSet.
+ // form so as to store it as a PresburgerSet.
IntMatrix activeRegionNorm = activeRegion.normalizeRows();
IntegerPolyhedron activeRegionRel = IntegerPolyhedron(
PresburgerSpace::getRelationSpace(0, numSymbols, 0, 0));
@@ -407,7 +404,7 @@ mlir::presburger::detail::computePolytopeGeneratingFunction(PolyhedronH poly) {
// 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.
+ // up to numVars - 1 dimensions.
std::vector<std::pair<PresburgerSet, std::vector<unsigned>>> chambers =
computeChamberDecomposition(activeRegions, vertices);
diff --git a/mlir/lib/Analysis/Presburger/IntegerRelation.cpp b/mlir/lib/Analysis/Presburger/IntegerRelation.cpp
index ab50f0dbec2cf..526933d119187 100644
--- a/mlir/lib/Analysis/Presburger/IntegerRelation.cpp
+++ b/mlir/lib/Analysis/Presburger/IntegerRelation.cpp
@@ -2516,10 +2516,10 @@ bool IntegerRelation::isFullDim() {
if (getNumEqualities() > 0)
return false;
- // If along the direction of any of the inequalities, the upper and lower
+ // If along the direction of any of the inequalities, the set is flat,
// optima are the same, then the region is not full-dimensional.
Simplex simplex(*this);
- return llvm::none_of(llvm::seq<int>(0, getNumInequalities()), [&](int i) {
+ return llvm::none_of(llvm::seq<int>(getNumInequalities()), [&](int i) {
return simplex.isFlatAlong(getInequality(i));
});
}
diff --git a/mlir/lib/Analysis/Presburger/PresburgerRelation.cpp b/mlir/lib/Analysis/Presburger/PresburgerRelation.cpp
index ee56ce696e86e..3af6baae0e700 100644
--- a/mlir/lib/Analysis/Presburger/PresburgerRelation.cpp
+++ b/mlir/lib/Analysis/Presburger/PresburgerRelation.cpp
@@ -1042,10 +1042,9 @@ PresburgerRelation PresburgerRelation::simplify() const {
}
bool PresburgerRelation::isFullDim() const {
- bool fullDim = llvm::any_of(getAllDisjuncts(), [&](IntegerRelation disjunct) {
+ return llvm::any_of(getAllDisjuncts(), [&](IntegerRelation disjunct) {
return disjunct.isFullDim();
});
- return fullDim;
}
void PresburgerRelation::print(raw_ostream &os) const {
diff --git a/mlir/lib/Analysis/Presburger/Simplex.cpp b/mlir/lib/Analysis/Presburger/Simplex.cpp
index 592d96ee766fb..1969cce93ad2e 100644
--- a/mlir/lib/Analysis/Presburger/Simplex.cpp
+++ b/mlir/lib/Analysis/Presburger/Simplex.cpp
@@ -2114,7 +2114,6 @@ bool Simplex::isFlatAlong(ArrayRef<MPInt> coeffs) {
if (!downOpt.isBounded())
return false;
- // Check if the upper and lower optima are equal.
return *upOpt == *downOpt;
}
>From 647a85461ac253c113a60d93341738d70f3c3f46 Mon Sep 17 00:00:00 2001
From: Abhinav271828 <abhinav.m at research.iiit.ac.in>
Date: Thu, 25 Jan 2024 12:40:59 +0530
Subject: [PATCH 30/41] Add negateMatrix and scaleRow
---
mlir/include/mlir/Analysis/Presburger/Matrix.h | 6 ++++++
mlir/lib/Analysis/Presburger/Barvinok.cpp | 11 +++--------
mlir/lib/Analysis/Presburger/Matrix.cpp | 12 ++++++++++++
3 files changed, 21 insertions(+), 8 deletions(-)
diff --git a/mlir/include/mlir/Analysis/Presburger/Matrix.h b/mlir/include/mlir/Analysis/Presburger/Matrix.h
index fe260568a22d8..83f642fba8d56 100644
--- a/mlir/include/mlir/Analysis/Presburger/Matrix.h
+++ b/mlir/include/mlir/Analysis/Presburger/Matrix.h
@@ -143,6 +143,9 @@ class Matrix {
/// Add `scale` multiples of the rowVec row to the specified row.
void addToRow(unsigned row, ArrayRef<T> rowVec, const T &scale);
+ /// Multiply the specified row by a factor of `scale`.
+ void scaleRow(unsigned row, const T &scale);
+
/// Add `scale` multiples of the source column to the target column.
void addToColumn(unsigned sourceColumn, unsigned targetColumn,
const T &scale);
@@ -157,6 +160,9 @@ class Matrix {
/// Negate the specified row.
void negateRow(unsigned row);
+ /// Negate the entire matrix.
+ void negateMatrix();
+
/// The given vector is interpreted as a row vector v. Post-multiply v with
/// this matrix, say M, and return vM.
SmallVector<T, 8> preMultiplyWithRow(ArrayRef<T> rowVec) const;
diff --git a/mlir/lib/Analysis/Presburger/Barvinok.cpp b/mlir/lib/Analysis/Presburger/Barvinok.cpp
index e289f3884f593..7419debf8afb1 100644
--- a/mlir/lib/Analysis/Presburger/Barvinok.cpp
+++ b/mlir/lib/Analysis/Presburger/Barvinok.cpp
@@ -204,11 +204,8 @@ mlir::presburger::detail::solveParametricEquations(FracMatrix equations) {
}
// Rescale diagonal elements to 1.
- for (unsigned i = 0; i < d; ++i) {
- Fraction diagElement = equations(i, i);
- for (unsigned j = 0; j < numCols; ++j)
- equations(i, j) = equations(i, j) / diagElement;
- }
+ for (unsigned i = 0; i < d; ++i)
+ equations.scaleRow(i, 1 / equations(i, i));
// Now we have reduced the equations to the form
// x_i + b_1' m_1 + ... + b_p' m_p + c' = 0
@@ -221,9 +218,7 @@ mlir::presburger::detail::solveParametricEquations(FracMatrix equations) {
ParamPoint vertex =
equations.getSubMatrix(/*fromRow=*/0, /*toRow=*/d - 1,
/*fromColumn=*/d, /*toColumn=*/numCols - 1);
- for (unsigned i = 0; i < d; ++i)
- vertex.negateRow(i);
-
+ vertex.negateMatrix();
return vertex;
}
diff --git a/mlir/lib/Analysis/Presburger/Matrix.cpp b/mlir/lib/Analysis/Presburger/Matrix.cpp
index 6d734c232f182..6f9697f1b4d28 100644
--- a/mlir/lib/Analysis/Presburger/Matrix.cpp
+++ b/mlir/lib/Analysis/Presburger/Matrix.cpp
@@ -295,6 +295,12 @@ void Matrix<T>::addToRow(unsigned row, ArrayRef<T> rowVec, const T &scale) {
at(row, col) += scale * rowVec[col];
}
+template <typename T>
+void Matrix<T>::scaleRow(unsigned row, const T &scale) {
+ for (unsigned col = 0; col < nColumns; ++col)
+ at(row, col) *= scale;
+}
+
template <typename T>
void Matrix<T>::addToColumn(unsigned sourceColumn, unsigned targetColumn,
const T &scale) {
@@ -316,6 +322,12 @@ void Matrix<T>::negateRow(unsigned row) {
at(row, column) = -at(row, column);
}
+template <typename T>
+void Matrix<T>::negateMatrix() {
+ for (unsigned row = 0; row < nRows; ++row)
+ negateRow(row);
+}
+
template <typename T>
SmallVector<T, 8> Matrix<T>::preMultiplyWithRow(ArrayRef<T> rowVec) const {
assert(rowVec.size() == getNumRows() && "Invalid row vector dimension!");
>From 06193a6c98c8b070ba4d794c6378ae4bee6d2226 Mon Sep 17 00:00:00 2001
From: Abhinav271828 <abhinav.m at research.iiit.ac.in>
Date: Thu, 25 Jan 2024 20:21:01 +0530
Subject: [PATCH 31/41] Update doc
---
mlir/lib/Analysis/Presburger/IntegerRelation.cpp | 4 ++--
1 file changed, 2 insertions(+), 2 deletions(-)
diff --git a/mlir/lib/Analysis/Presburger/IntegerRelation.cpp b/mlir/lib/Analysis/Presburger/IntegerRelation.cpp
index 526933d119187..2ac271e2e0553 100644
--- a/mlir/lib/Analysis/Presburger/IntegerRelation.cpp
+++ b/mlir/lib/Analysis/Presburger/IntegerRelation.cpp
@@ -2516,8 +2516,8 @@ bool IntegerRelation::isFullDim() {
if (getNumEqualities() > 0)
return false;
- // If along the direction of any of the inequalities, the set is flat,
- // optima are the same, then the region is not full-dimensional.
+ // The polytope is full-dimensional iff it is not flat along any of the
+ // inequality directions.
Simplex simplex(*this);
return llvm::none_of(llvm::seq<int>(getNumInequalities()), [&](int i) {
return simplex.isFlatAlong(getInequality(i));
>From b7f9ba3bf014e20aa78468e635cd056062c53595 Mon Sep 17 00:00:00 2001
From: Abhinav271828 <abhinav.m at research.iiit.ac.in>
Date: Thu, 25 Jan 2024 20:47:44 +0530
Subject: [PATCH 32/41] Refactor
---
.../mlir/Analysis/Presburger/Barvinok.h | 7 +-
mlir/lib/Analysis/Presburger/Barvinok.cpp | 164 +++++++++---------
.../Analysis/Presburger/BarvinokTest.cpp | 21 +--
3 files changed, 88 insertions(+), 104 deletions(-)
diff --git a/mlir/include/mlir/Analysis/Presburger/Barvinok.h b/mlir/include/mlir/Analysis/Presburger/Barvinok.h
index 389be6914f3cc..0974378e09b66 100644
--- a/mlir/include/mlir/Analysis/Presburger/Barvinok.h
+++ b/mlir/include/mlir/Analysis/Presburger/Barvinok.h
@@ -105,9 +105,10 @@ std::optional<ParamPoint> solveParametricEquations(FracMatrix equations);
/// The ith relation corresponds to the activity region of the ith vertex set.
/// Note that here, by disjoint, we mean that the intersection is not
/// full-dimensional.
-std::vector<std::pair<PresburgerSet, std::vector<unsigned>>>
-computeChamberDecomposition(ArrayRef<PresburgerSet> activeRegions,
- ArrayRef<ParamPoint> vertices);
+std::vector<std::pair<PresburgerSet, GeneratingFunction>>
+computeChamberDecomposition(
+ unsigned numSymbols, ArrayRef<std::pair<PresburgerSet, GeneratingFunction>>
+ regionsAndGeneratingFunctions);
/// Compute the generating function corresponding to a polytope.
/// All tangent cones of the polytope must be unimodular.
diff --git a/mlir/lib/Analysis/Presburger/Barvinok.cpp b/mlir/lib/Analysis/Presburger/Barvinok.cpp
index 7419debf8afb1..0034753001d53 100644
--- a/mlir/lib/Analysis/Presburger/Barvinok.cpp
+++ b/mlir/lib/Analysis/Presburger/Barvinok.cpp
@@ -224,62 +224,70 @@ mlir::presburger::detail::solveParametricEquations(FracMatrix equations) {
/// 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 separating the chambers where it is not.
+/// We maintain a list of pairwise disjoint chambers and their generating
+/// functions. We iterate over the list of regions, each time adding the
+/// generating function to the chambers where the corresponding vertex is
+/// active and separating the chambers where it is not.
///
/// Given the region each vertex is active in, for each subset of vertices,
/// the region that precisely this subset is in, is the intersection of the
/// regions that these are active in, intersected with the complements of the
-/// remaining regions.
-std::vector<std::pair<PresburgerSet, std::vector<unsigned>>>
+/// remaining regions. The generating function for this region is the sum of
+/// generating functions for the vertices' tangent cones.
+std::vector<std::pair<PresburgerSet, GeneratingFunction>>
mlir::presburger::detail::computeChamberDecomposition(
- ArrayRef<PresburgerSet> activeRegions, ArrayRef<ParamPoint> vertices) {
- // We maintain a list of regions and their associated vertex sets,
- // initialized with the universe and the empty list of vertices.
- std::vector<std::pair<PresburgerSet, std::vector<unsigned>>> chambers = {
- {PresburgerSet::getUniverse(activeRegions[0].getSpace()),
- std::vector<unsigned>({})}};
- // 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]}.
-
- // We iterate over the vertex set.
- // For each vertex v_j and its activity region R_j,
+ unsigned numSymbols, ArrayRef<std::pair<PresburgerSet, GeneratingFunction>>
+ regionsAndGeneratingFunctions) {
+ assert(regionsAndGeneratingFunctions.size() > 0 &&
+ "there must be at least one chamber!");
+ // We maintain a list of regions and their associated generating function
+ // initialized with the universe and the empty generating function.
+ std::vector<std::pair<PresburgerSet, GeneratingFunction>> chambers = {
+ {PresburgerSet::getUniverse(
+ PresburgerSpace::getRelationSpace(0, numSymbols, 0, 0)),
+ GeneratingFunction(numSymbols, {}, {}, {})}};
+
+ // We iterate over the region list.
+ // For each activity region R_j (corresponding to a vertex v_j, whose
+ // generating function is gf_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.
+ // 1. the intersection R_i \cap R_j, where the generating function is
+ // gf_i + gf_j.
+ // 2. the difference R_i - R_j, where v_j is inactive and the generating
+ // function is gf_i.
// At each step, we define a new chamber list after considering vertex v_j,
- // replacing and appending chambers as discussed above.
+ // and its generating function, replacing and appending chambers as
+ // discussed above.
// The loop has the invariant that the union over all the chambers gives the
// universe at every step (modulo lower-dimensional spaces).
- for (unsigned j = 0, e = vertices.size(); j < e; j++) {
- std::vector<std::pair<PresburgerSet, std::vector<unsigned>>> newChambers;
+ for (const auto &[region, generatingFunction] :
+ regionsAndGeneratingFunctions) {
+ std::vector<std::pair<PresburgerSet, GeneratingFunction>> newChambers;
- for (auto [currentRegion, currentVertices] : chambers) {
+ for (auto [currentRegion, currentGeneratingFunction] : chambers) {
// First, we check if the intersection of R_j and R_i.
- PresburgerSet intersection = currentRegion.intersect(activeRegions[j]);
+ PresburgerSet intersection = currentRegion.intersect(region);
// If the intersection is not full-dimensional, we do not modify
// the chamber list.
if (!intersection.isFullDim()) {
- newChambers.emplace_back(currentRegion, currentVertices);
+ newChambers.emplace_back(currentRegion, currentGeneratingFunction);
continue;
}
// If it is, we add the intersection as a chamber where this vertex is
- // active. We also add the difference if it is full-dimensional.
- PresburgerSet difference = currentRegion.subtract(activeRegions[j]);
- if (difference.isFullDim())
- newChambers.emplace_back(difference, currentVertices);
+ // active.
+ newChambers.emplace_back(intersection,
+ currentGeneratingFunction + generatingFunction);
- currentVertices.push_back(j);
- newChambers.emplace_back(intersection, currentVertices);
+ // We also add the difference if it is full-dimensional.
+ PresburgerSet difference = currentRegion.subtract(region);
+ if (difference.isFullDim())
+ newChambers.emplace_back(difference, currentGeneratingFunction);
}
-
chambers = std::move(newChambers);
}
@@ -307,17 +315,15 @@ mlir::presburger::detail::computePolytopeGeneratingFunction(PolyhedronH poly) {
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<PresburgerSet, GeneratingFunction>> gf({});
-
// 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;
+ // the vertices corresponding to them,
+ // the active regions of the vertices, and the generating functions of
+ // the tangent cones of the vertices, in order.
+ // std::vector<IntMatrix> subsets;
std::vector<ParamPoint> vertices;
- std::vector<PresburgerSet> activeRegions;
+ std::vector<std::pair<PresburgerSet, GeneratingFunction>>
+ regionsAndGeneratingFunctions;
FracMatrix a2(numIneqs - numVars, numVars);
FracMatrix b2c2(numIneqs - numVars, numSymbols + 1);
@@ -360,7 +366,7 @@ mlir::presburger::detail::computePolytopeGeneratingFunction(PolyhedronH poly) {
continue;
// If this subset corresponds to a vertex, store it.
vertices.push_back(*vertex);
- subsets.push_back(subset);
+ // subsets.push_back(subset);
// Let the current vertex be [X | y], where
// X represents the coefficients of the parameters and
@@ -392,7 +398,38 @@ mlir::presburger::detail::computePolytopeGeneratingFunction(PolyhedronH poly) {
PresburgerSpace::getRelationSpace(0, numSymbols, 0, 0));
for (unsigned i = 0, e = activeRegion.getNumRows(); i < e; ++i)
activeRegionRel.addInequality(activeRegionNorm.getRow(i));
- activeRegions.push_back(PresburgerSet(activeRegionRel));
+ // TODO check for simplicial
+
+ // Now, we compute the generating function at this vertex.
+ // We collect the inequalities corresponding to each vertex to compute
+ // the tangent cone at that vertex.
+ SmallVector<MPInt> ineq(numVars + 1);
+
+ // We only need the coefficients of the variables (NOT the parameters)
+ // as the generating function only depends on these.
+ // We translate the cones to be pointed at the origin by making the
+ // constant terms zero.
+ ConeH tangentCone = defineHRep(numVars);
+ for (unsigned j = 0; j < numVars; ++j) {
+ for (unsigned k = 0; k < numVars; ++k)
+ ineq[k] = subset(j, k);
+ tangentCone.addInequality(ineq);
+ }
+ // We assume that the tangent cone is unimodular, so there is no need
+ // to decompose it.
+ // In the general case, the unimodular decomposition may have several
+ // cones.
+ GeneratingFunction vertexGf(numSymbols, {}, {}, {});
+ SmallVector<std::pair<int, ConeH>, 4> unimodCones = {{1, tangentCone}};
+ for (std::pair<int, ConeH> signedCone : unimodCones) {
+ auto [sign, cone] = signedCone;
+ vertexGf = vertexGf +
+ computeUnimodularConeGeneratingFunction(*vertex, sign, cone);
+ }
+ // We store the vertex we computed with the generating function of its
+ // tangent cones.
+ regionsAndGeneratingFunctions.emplace_back(PresburgerSet(activeRegionRel),
+ vertexGf);
}
// Now, we use Clauss-Loechner decomposition to identify regions in parameter
@@ -400,46 +437,9 @@ mlir::presburger::detail::computePolytopeGeneratingFunction(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.
- std::vector<std::pair<PresburgerSet, std::vector<unsigned>>> chambers =
- computeChamberDecomposition(activeRegions, vertices);
-
- // Now, we compute the generating function. For each chamber, we iterate over
- // the vertices active in it, and compute the generating function for each
- // 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<PresburgerSet, std::vector<unsigned>> &chamber :
- chambers) {
- auto [currentRegion, currentVertices] = chamber;
- GeneratingFunction chamberGf(numSymbols, {}, {}, {});
- for (unsigned i : currentVertices) {
- // We collect the inequalities corresponding to each vertex to compute
- // the tangent cone at that vertex.
- // We only need the coefficients of the variables (NOT the parameters)
- // as the generating function only depends on these.
- // We translate the cones to be pointed at the origin by making the
- // constant terms zero.
- ConeH tangentCone = defineHRep(numVars);
- for (unsigned j = 0; j < numVars; ++j) {
- for (unsigned k = 0; k < numVars; ++k)
- ineq[k] = subsets[i](j, k);
- tangentCone.addInequality(ineq);
- }
- // We assume that the tangent cone is unimodular, so there is no need
- // to decompose it.
- // In the general case, the unimodular decomposition may have several
- // cones.
- SmallVector<std::pair<int, ConeH>, 4> unimodCones = {
- std::make_pair(1, tangentCone)};
- for (std::pair<int, ConeH> signedCone : unimodCones) {
- auto [sign, cone] = signedCone;
- chamberGf = chamberGf + computeUnimodularConeGeneratingFunction(
- vertices[i], sign, cone);
- }
- }
- gf.emplace_back(currentRegion, chamberGf);
- }
- return gf;
+ // In each chamber, we sum up the generating functions of the active vertices
+ // to find the generating function of the polytope.
+ return computeChamberDecomposition(numSymbols, regionsAndGeneratingFunctions);
}
/// We use an iterative procedure to find a vector not orthogonal
diff --git a/mlir/unittests/Analysis/Presburger/BarvinokTest.cpp b/mlir/unittests/Analysis/Presburger/BarvinokTest.cpp
index 2ab39d0e10ea5..40219c3ec7a69 100644
--- a/mlir/unittests/Analysis/Presburger/BarvinokTest.cpp
+++ b/mlir/unittests/Analysis/Presburger/BarvinokTest.cpp
@@ -286,7 +286,7 @@ TEST(BarvinokTest, computeNumTermsPolytope) {
// There is only one chamber: p ≥ 0
EXPECT_EQ(count.size(), 2u);
- gf = count[1].second;
+ gf = count[0].second;
EXPECT_EQ_REPR_GENERATINGFUNCTION(
gf, GeneratingFunction(
1, {1, 1, 1},
@@ -295,23 +295,6 @@ TEST(BarvinokTest, computeNumTermsPolytope) {
makeFracMatrix(2, 2, {{0, 0}, {0, 0}})},
{{{-1, 1}, {-1, 0}}, {{1, -1}, {0, -1}}, {{1, 0}, {0, 1}}}));
- ineqs = makeIntMatrix(6, 5,
- {{1, 0, 0, 1, 0},
- {0, 1, 0, 1, 0},
- {0, 0, 1, 1, 0},
- {-1, 0, 0, 1, 0},
- {0, -1, 0, 1, 0},
- {0, 0, -1, 1, 0}});
- poly = defineHRep(3, 1);
- for (unsigned i = 0; i < 6; i++)
- poly.addInequality(ineqs.getRow(i));
-
- count = computePolytopeGeneratingFunction(poly);
-
- EXPECT_EQ(count.size(), 2u);
-
- gf = count[0].second;
-
// Cartesian product of a cube with side M and a right triangle with side N.
ineqs = makeIntMatrix(9, 8,
{{1, 0, 0, 0, 0, 1, 0, 0},
@@ -331,6 +314,6 @@ TEST(BarvinokTest, computeNumTermsPolytope) {
EXPECT_EQ(count.size(), 2u);
- gf = count[1].second;
+ gf = count[0].second;
EXPECT_EQ(gf.getNumerators().size(), 24u);
}
>From ec1023e083b3310bb9b0f8a4b22bda06ed2cba56 Mon Sep 17 00:00:00 2001
From: Abhinav271828 <abhinav.m at research.iiit.ac.in>
Date: Thu, 25 Jan 2024 20:56:18 +0530
Subject: [PATCH 33/41] Fix nonsimplicial case
---
mlir/lib/Analysis/Presburger/Barvinok.cpp | 19 +++++++++++++++++--
1 file changed, 17 insertions(+), 2 deletions(-)
diff --git a/mlir/lib/Analysis/Presburger/Barvinok.cpp b/mlir/lib/Analysis/Presburger/Barvinok.cpp
index 0034753001d53..97124f5600ea5 100644
--- a/mlir/lib/Analysis/Presburger/Barvinok.cpp
+++ b/mlir/lib/Analysis/Presburger/Barvinok.cpp
@@ -391,6 +391,22 @@ mlir::presburger::detail::computePolytopeGeneratingFunction(PolyhedronH poly) {
activeRegion.addToRow(i, b2c2.getRow(i), 1);
}
+ // If any row of activeRegion is all-zero, that means that the
+ // corresponding inequality in `remaining` is *also* satisfied by the
+ // vertex for any values of the parameters. Thus we remove it from
+ // activeRegion and add that inequality to `subset`.
+ // Note that if the row is not all-zero, it is still possible for the
+ // inequality to be satisfied by the corresponding vertex *in some subset of
+ // the parameter space*. However, since this subset is defined by an
+ // equality, it is not full-dimensional and we can therefore ignore it.
+ for (int i = activeRegion.getNumRows() - 1; i >= 0; --i) {
+ if (llvm::any_of(activeRegion.getRow(i),
+ [&](Fraction f) { return f != 0; }))
+ continue;
+ subset.appendExtraRow(remainder.getRow(i));
+ activeRegion.removeRow(i);
+ }
+
// We convert the representation of the active region to an integers-only
// form so as to store it as a PresburgerSet.
IntMatrix activeRegionNorm = activeRegion.normalizeRows();
@@ -398,7 +414,6 @@ mlir::presburger::detail::computePolytopeGeneratingFunction(PolyhedronH poly) {
PresburgerSpace::getRelationSpace(0, numSymbols, 0, 0));
for (unsigned i = 0, e = activeRegion.getNumRows(); i < e; ++i)
activeRegionRel.addInequality(activeRegionNorm.getRow(i));
- // TODO check for simplicial
// Now, we compute the generating function at this vertex.
// We collect the inequalities corresponding to each vertex to compute
@@ -410,7 +425,7 @@ mlir::presburger::detail::computePolytopeGeneratingFunction(PolyhedronH poly) {
// We translate the cones to be pointed at the origin by making the
// constant terms zero.
ConeH tangentCone = defineHRep(numVars);
- for (unsigned j = 0; j < numVars; ++j) {
+ for (unsigned j = 0, e = subset.getNumRows(); j < e; ++j) {
for (unsigned k = 0; k < numVars; ++k)
ineq[k] = subset(j, k);
tangentCone.addInequality(ineq);
>From 6ffcaca6d23cf8a24e8d1db560f7334cc0234477 Mon Sep 17 00:00:00 2001
From: Abhinav271828 <abhinav.m at research.iiit.ac.in>
Date: Fri, 26 Jan 2024 19:22:26 +0530
Subject: [PATCH 34/41] Fixes
---
.../mlir/Analysis/Presburger/Barvinok.h | 2 +-
.../Analysis/Presburger/IntegerRelation.h | 9 +++
mlir/include/mlir/Analysis/Presburger/Utils.h | 2 +
mlir/lib/Analysis/Presburger/Barvinok.cpp | 67 +++++++++----------
mlir/lib/Analysis/Presburger/Utils.cpp | 4 ++
5 files changed, 47 insertions(+), 37 deletions(-)
diff --git a/mlir/include/mlir/Analysis/Presburger/Barvinok.h b/mlir/include/mlir/Analysis/Presburger/Barvinok.h
index 0974378e09b66..a80aa9bc2f9b2 100644
--- a/mlir/include/mlir/Analysis/Presburger/Barvinok.h
+++ b/mlir/include/mlir/Analysis/Presburger/Barvinok.h
@@ -113,7 +113,7 @@ computeChamberDecomposition(
/// Compute the generating function corresponding to a polytope.
/// All tangent cones of the polytope must be unimodular.
std::vector<std::pair<PresburgerSet, GeneratingFunction>>
-computePolytopeGeneratingFunction(PolyhedronH poly);
+computePolytopeGeneratingFunction(const 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/include/mlir/Analysis/Presburger/IntegerRelation.h b/mlir/include/mlir/Analysis/Presburger/IntegerRelation.h
index 4c4576bf6d0d3..b4159b490381e 100644
--- a/mlir/include/mlir/Analysis/Presburger/IntegerRelation.h
+++ b/mlir/include/mlir/Analysis/Presburger/IntegerRelation.h
@@ -882,6 +882,15 @@ class IntegerPolyhedron : public IntegerRelation {
/*numReservedEqualities=*/0,
/*numReservedCols=*/space.getNumVars() + 1, space) {}
+ /// Constructs a relation with the specified number of dimensions and symbols
+ /// and adds the given inequalities.
+ explicit IntegerPolyhedron(const PresburgerSpace &space,
+ IntMatrix inequalities)
+ : IntegerPolyhedron(space) {
+ for (unsigned i = 0, e = inequalities.getNumRows(); i < e; i++)
+ addInequality(inequalities.getRow(i));
+ }
+
/// Construct a set from an IntegerRelation. The relation should have
/// no domain vars.
explicit IntegerPolyhedron(const IntegerRelation &rel)
diff --git a/mlir/include/mlir/Analysis/Presburger/Utils.h b/mlir/include/mlir/Analysis/Presburger/Utils.h
index e6d29e4ef6d06..38262a65f9754 100644
--- a/mlir/include/mlir/Analysis/Presburger/Utils.h
+++ b/mlir/include/mlir/Analysis/Presburger/Utils.h
@@ -286,6 +286,8 @@ Fraction dotProduct(ArrayRef<Fraction> a, ArrayRef<Fraction> b);
std::vector<Fraction> multiplyPolynomials(ArrayRef<Fraction> a,
ArrayRef<Fraction> b);
+bool isRangeZero(ArrayRef<Fraction> arr);
+
} // namespace presburger
} // namespace mlir
diff --git a/mlir/lib/Analysis/Presburger/Barvinok.cpp b/mlir/lib/Analysis/Presburger/Barvinok.cpp
index 97124f5600ea5..ecb88064ed1b7 100644
--- a/mlir/lib/Analysis/Presburger/Barvinok.cpp
+++ b/mlir/lib/Analysis/Presburger/Barvinok.cpp
@@ -238,25 +238,23 @@ std::vector<std::pair<PresburgerSet, GeneratingFunction>>
mlir::presburger::detail::computeChamberDecomposition(
unsigned numSymbols, ArrayRef<std::pair<PresburgerSet, GeneratingFunction>>
regionsAndGeneratingFunctions) {
- assert(regionsAndGeneratingFunctions.size() > 0 &&
+ assert(!regionsAndGeneratingFunctions.empty() &&
"there must be at least one chamber!");
// We maintain a list of regions and their associated generating function
// initialized with the universe and the empty generating function.
std::vector<std::pair<PresburgerSet, GeneratingFunction>> chambers = {
- {PresburgerSet::getUniverse(
- PresburgerSpace::getRelationSpace(0, numSymbols, 0, 0)),
+ {PresburgerSet::getUniverse(PresburgerSpace::getSetSpace(numSymbols)),
GeneratingFunction(numSymbols, {}, {}, {})}};
// We iterate over the region list.
// For each activity region R_j (corresponding to a vertex v_j, whose
- // generating function is gf_j),
- // we examine all the current chambers R_i.
+ // generating function is gf_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 the generating function is
- // gf_i + gf_j.
- // 2. the difference R_i - R_j, where v_j is inactive and the generating
- // function is gf_i.
+ // 1. the intersection R_i \cap R_j (if it is full-dimensional), where the
+ // generating function is gf_i + gf_j.
+ // 2. the difference R_i - R_j (if it is full-dimensional), where v_j is
+ // inactive and the generating function is gf_i.
// At each step, we define a new chamber list after considering vertex v_j,
// and its generating function, replacing and appending chambers as
@@ -267,8 +265,7 @@ mlir::presburger::detail::computeChamberDecomposition(
regionsAndGeneratingFunctions) {
std::vector<std::pair<PresburgerSet, GeneratingFunction>> newChambers;
- for (auto [currentRegion, currentGeneratingFunction] : chambers) {
- // First, we check if the intersection of R_j and R_i.
+ for (const auto &[currentRegion, currentGeneratingFunction] : chambers) {
PresburgerSet intersection = currentRegion.intersect(region);
// If the intersection is not full-dimensional, we do not modify
@@ -294,14 +291,18 @@ mlir::presburger::detail::computeChamberDecomposition(
return chambers;
}
-/// For a polytope expressed as a set of inequalities, compute the generating
-/// function corresponding to the number of lattice points present. This
+/// For a polytope expressed as a set of n inequalities, compute the generating
+/// function corresponding to the lattice points included in the polytope. This
/// algorithm has three main steps:
/// 1. Enumerate the vertices, by iterating over subsets of inequalities and
-/// checking for satisfiability.
+/// checking for satisfiability. For each d-subset of inequalities (where d
+/// is the number of variables), we solve to obtain the vertex in terms of
+/// the parameters, and then check for the region in parameter space where
+/// this vertex satisfies the remaining (n - d) inequalities.
/// 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.
+/// function corresponding to it. The generating function depends on the
+/// parametric expression of the vertex and the (non-parametric) generators
+/// of the tangent cone.
/// 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.
@@ -310,24 +311,19 @@ mlir::presburger::detail::computeChamberDecomposition(
/// polytopes using Barvinok's rational functions." Algorithmica 48 (2007):
/// 37-66.
std::vector<std::pair<PresburgerSet, GeneratingFunction>>
-mlir::presburger::detail::computePolytopeGeneratingFunction(PolyhedronH poly) {
+mlir::presburger::detail::computePolytopeGeneratingFunction(
+ const PolyhedronH &poly) {
unsigned numVars = poly.getNumRangeVars();
unsigned numSymbols = poly.getNumSymbolVars();
unsigned numIneqs = poly.getNumInequalities();
- // These vectors store lists of
- // subsets of inequalities,
- // the vertices corresponding to them,
- // the active regions of the vertices, and the generating functions of
- // the tangent cones of the vertices, in order.
- // std::vector<IntMatrix> subsets;
+ // We store a list of the computed vertices.
std::vector<ParamPoint> vertices;
+ // For each vertex, we store the corresponding active region and the
+ // generating functions of the tangent cone, in order.
std::vector<std::pair<PresburgerSet, GeneratingFunction>>
regionsAndGeneratingFunctions;
- FracMatrix a2(numIneqs - numVars, numVars);
- FracMatrix b2c2(numIneqs - numVars, numSymbols + 1);
-
// We iterate over all subsets of inequalities with cardinality numVars,
// using bitsets of numIneqs bits to enumerate.
// For a given set of numIneqs bits, we consider a subset which contains
@@ -352,6 +348,8 @@ mlir::presburger::detail::computePolytopeGeneratingFunction(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.
+ FracMatrix a2(numIneqs - numVars, numVars);
+ FracMatrix b2c2(numIneqs - numVars, numSymbols + 1);
a2 = FracMatrix(
remainder.getSubMatrix(0, numIneqs - numVars - 1, 0, numVars - 1));
b2c2 = FracMatrix(remainder.getSubMatrix(0, numIneqs - numVars - 1, numVars,
@@ -366,7 +364,6 @@ mlir::presburger::detail::computePolytopeGeneratingFunction(PolyhedronH poly) {
continue;
// If this subset corresponds to a vertex, store it.
vertices.push_back(*vertex);
- // subsets.push_back(subset);
// Let the current vertex be [X | y], where
// X represents the coefficients of the parameters and
@@ -398,10 +395,10 @@ mlir::presburger::detail::computePolytopeGeneratingFunction(PolyhedronH poly) {
// Note that if the row is not all-zero, it is still possible for the
// inequality to be satisfied by the corresponding vertex *in some subset of
// the parameter space*. However, since this subset is defined by an
- // equality, it is not full-dimensional and we can therefore ignore it.
+ // equality, it is not full-dimensional; here and in the chamber
+ // decomposition, we ignore full-dimensional chambers.
for (int i = activeRegion.getNumRows() - 1; i >= 0; --i) {
- if (llvm::any_of(activeRegion.getRow(i),
- [&](Fraction f) { return f != 0; }))
+ if (!isRangeZero(activeRegion.getRow(i)))
continue;
subset.appendExtraRow(remainder.getRow(i));
activeRegion.removeRow(i);
@@ -409,11 +406,9 @@ mlir::presburger::detail::computePolytopeGeneratingFunction(PolyhedronH poly) {
// We convert the representation of the active region to an integers-only
// form so as to store it as a PresburgerSet.
- IntMatrix activeRegionNorm = activeRegion.normalizeRows();
- IntegerPolyhedron activeRegionRel = IntegerPolyhedron(
- PresburgerSpace::getRelationSpace(0, numSymbols, 0, 0));
- for (unsigned i = 0, e = activeRegion.getNumRows(); i < e; ++i)
- activeRegionRel.addInequality(activeRegionNorm.getRow(i));
+ IntegerPolyhedron activeRegionRel(
+ PresburgerSpace::getRelationSpace(0, numSymbols, 0, 0),
+ activeRegion.normalizeRows());
// Now, we compute the generating function at this vertex.
// We collect the inequalities corresponding to each vertex to compute
@@ -450,7 +445,7 @@ mlir::presburger::detail::computePolytopeGeneratingFunction(PolyhedronH poly) {
// 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
+ // they may share "facets" or "edges", but their intersection can only have
// up to numVars - 1 dimensions.
// In each chamber, we sum up the generating functions of the active vertices
// to find the generating function of the polytope.
diff --git a/mlir/lib/Analysis/Presburger/Utils.cpp b/mlir/lib/Analysis/Presburger/Utils.cpp
index a8d860885ef10..f89373d8ba506 100644
--- a/mlir/lib/Analysis/Presburger/Utils.cpp
+++ b/mlir/lib/Analysis/Presburger/Utils.cpp
@@ -564,4 +564,8 @@ std::vector<Fraction> presburger::multiplyPolynomials(ArrayRef<Fraction> a,
convolution.push_back(sum);
}
return convolution;
+}
+
+bool presburger::isRangeZero(ArrayRef<Fraction> arr) {
+ return llvm::all_of(arr, [&](Fraction f) { return f == 0; });
}
\ No newline at end of file
>From 6ce07a38a06a52d107155e9e88f0186bbb6b1db1 Mon Sep 17 00:00:00 2001
From: Abhinav271828 <abhinav.m at research.iiit.ac.in>
Date: Fri, 26 Jan 2024 22:19:20 +0530
Subject: [PATCH 35/41] Refactor
---
.../include/mlir/Analysis/Presburger/Matrix.h | 2 +-
mlir/lib/Analysis/Presburger/Barvinok.cpp | 27 +++----
mlir/lib/Analysis/Presburger/Matrix.cpp | 4 +-
.../Analysis/Presburger/BarvinokTest.cpp | 74 ++++++++++---------
4 files changed, 51 insertions(+), 56 deletions(-)
diff --git a/mlir/include/mlir/Analysis/Presburger/Matrix.h b/mlir/include/mlir/Analysis/Presburger/Matrix.h
index 83f642fba8d56..15b39a4fd927f 100644
--- a/mlir/include/mlir/Analysis/Presburger/Matrix.h
+++ b/mlir/include/mlir/Analysis/Presburger/Matrix.h
@@ -201,7 +201,7 @@ class Matrix {
/// 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);
+ std::pair<Matrix<T>, Matrix<T>> splitByBitset(std::vector<int> indicator);
/// Print the matrix.
void print(raw_ostream &os) const;
diff --git a/mlir/lib/Analysis/Presburger/Barvinok.cpp b/mlir/lib/Analysis/Presburger/Barvinok.cpp
index ecb88064ed1b7..fab8ac28af60e 100644
--- a/mlir/lib/Analysis/Presburger/Barvinok.cpp
+++ b/mlir/lib/Analysis/Presburger/Barvinok.cpp
@@ -325,22 +325,15 @@ mlir::presburger::detail::computePolytopeGeneratingFunction(
regionsAndGeneratingFunctions;
// We iterate over all subsets of inequalities with cardinality numVars,
- // 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(upperBound);
- indicator.to_ulong() <= upperBound;
- indicator = std::bitset<16>(indicator.to_ulong() - 1)) {
-
- if (indicator.count() != numVars)
- continue;
-
+ // using permutations of numVars 1's and (numIneqs - numVars) 0's.
+ // For a given permutation, we consider a subset which contains
+ // the i'th inequality if the i'th bit in the bitset is 1.
+ // We start with the permutation that takes the last numVars inequalities.
+ std::vector<int> indicator(numIneqs);
+ for (unsigned i = numIneqs - numVars; i < numIneqs; ++i)
+ indicator[i] = 1;
+
+ do {
// Collect the inequalities corresponding to the bits which are set
// and the remaining ones.
auto [subset, remainder] = poly.getInequalities().splitByBitset(indicator);
@@ -440,7 +433,7 @@ mlir::presburger::detail::computePolytopeGeneratingFunction(
// tangent cones.
regionsAndGeneratingFunctions.emplace_back(PresburgerSet(activeRegionRel),
vertexGf);
- }
+ } while (std::next_permutation(indicator.begin(), indicator.end()));
// Now, we use Clauss-Loechner decomposition to identify regions in parameter
// space where each vertex is active. These regions (chambers) have the
diff --git a/mlir/lib/Analysis/Presburger/Matrix.cpp b/mlir/lib/Analysis/Presburger/Matrix.cpp
index 6f9697f1b4d28..3b251a1ca2266 100644
--- a/mlir/lib/Analysis/Presburger/Matrix.cpp
+++ b/mlir/lib/Analysis/Presburger/Matrix.cpp
@@ -395,10 +395,10 @@ void Matrix<T>::print(raw_ostream &os) const {
/// 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>::splitByBitset(std::vector<int> indicator) {
Matrix<T> rowsForOne(0, nColumns), rowsForZero(0, nColumns);
for (unsigned i = 0; i < nRows; i++) {
- if (indicator.test(i))
+ if (indicator[i] == 1)
rowsForOne.appendExtraRow(getRow(i));
else
rowsForZero.appendExtraRow(getRow(i));
diff --git a/mlir/unittests/Analysis/Presburger/BarvinokTest.cpp b/mlir/unittests/Analysis/Presburger/BarvinokTest.cpp
index 40219c3ec7a69..6799dd37caf1f 100644
--- a/mlir/unittests/Analysis/Presburger/BarvinokTest.cpp
+++ b/mlir/unittests/Analysis/Presburger/BarvinokTest.cpp
@@ -51,38 +51,40 @@ TEST(BarvinokTest, getIndex) {
// (s.t. the cones are unimodular) and the generating functions
// are computed. We check that the results contain the correct
// matrices.
-TEST(BarvinokTest, unimodularConeGeneratingFunction) {
- ConeH cone = defineHRep(2);
- cone.addInequality({0, -1, 0});
- cone.addInequality({-1, -2, 0});
-
- ParamPoint vertex =
- makeFracMatrix(2, 3, {{2, 2, 0}, {-1, -Fraction(1, 2), 1}});
-
- GeneratingFunction gf =
- computeUnimodularConeGeneratingFunction(vertex, 1, cone);
-
- EXPECT_EQ_REPR_GENERATINGFUNCTION(
- gf, GeneratingFunction(
- 2, {1},
- {makeFracMatrix(3, 2, {{-1, 0}, {-Fraction(1, 2), 1}, {1, 2}})},
- {{{2, -1}, {-1, 0}}}));
-
- cone = defineHRep(3);
- cone.addInequality({7, 1, 6, 0});
- cone.addInequality({9, 1, 7, 0});
- cone.addInequality({8, -1, 1, 0});
-
- vertex = makeFracMatrix(3, 2, {{5, 2}, {6, 2}, {7, 1}});
-
- gf = computeUnimodularConeGeneratingFunction(vertex, 1, cone);
-
- EXPECT_EQ_REPR_GENERATINGFUNCTION(
- gf,
- GeneratingFunction(
- 1, {1}, {makeFracMatrix(2, 3, {{-83, -100, -41}, {-22, -27, -15}})},
- {{{8, 47, -17}, {-7, -41, 15}, {1, 5, -2}}}));
-}
+// TEST(BarvinokTest, unimodularConeGeneratingFunction) {
+// ConeH cone = defineHRep(2);
+// cone.addInequality({0, -1, 0});
+// cone.addInequality({-1, -2, 0});
+//
+// ParamPoint vertex =
+// makeFracMatrix(2, 3, {{2, 2, 0}, {-1, -Fraction(1, 2), 1}});
+//
+// GeneratingFunction gf =
+// computeUnimodularConeGeneratingFunction(vertex, 1, cone);
+//
+// EXPECT_EQ_REPR_GENERATINGFUNCTION(
+// gf, GeneratingFunction(
+// 2, {1},
+// {makeFracMatrix(3, 2, {{-1, 0}, {-Fraction(1, 2), 1}, {1,
+// 2}})},
+// {{{2, -1}, {-1, 0}}}));
+//
+// cone = defineHRep(3);
+// cone.addInequality({7, 1, 6, 0});
+// cone.addInequality({9, 1, 7, 0});
+// cone.addInequality({8, -1, 1, 0});
+//
+// vertex = makeFracMatrix(3, 2, {{5, 2}, {6, 2}, {7, 1}});
+//
+// gf = computeUnimodularConeGeneratingFunction(vertex, 1, cone);
+//
+// EXPECT_EQ_REPR_GENERATINGFUNCTION(
+// gf,
+// GeneratingFunction(
+// 1, {1}, {makeFracMatrix(2, 3, {{-83, -100, -41}, {-22, -27,
+// -15}})},
+// {{{8, 47, -17}, {-7, -41, 15}, {1, 5, -2}}}));
+// }
// The following vectors are randomly generated.
// We then check that the output of the function has non-zero
@@ -268,12 +270,12 @@ TEST(BarvinokTest, computeNumTermsPolytope) {
makeFracMatrix(1, 3, {{0, 0, 1}}),
makeFracMatrix(1, 3, {{0, 0, 0}})},
{{{-1, 0, 0}, {0, -1, 0}, {0, 0, -1}},
- {{1, 0, 0}, {0, -1, 0}, {0, 0, -1}},
- {{0, 1, 0}, {-1, 0, 0}, {0, 0, -1}},
- {{1, 0, 0}, {0, 1, 0}, {0, 0, -1}},
{{0, 0, 1}, {-1, 0, 0}, {0, -1, 0}},
- {{1, 0, 0}, {0, 0, 1}, {0, -1, 0}},
+ {{0, 1, 0}, {-1, 0, 0}, {0, 0, -1}},
{{0, 1, 0}, {0, 0, 1}, {-1, 0, 0}},
+ {{1, 0, 0}, {0, -1, 0}, {0, 0, -1}},
+ {{1, 0, 0}, {0, 0, 1}, {0, -1, 0}},
+ {{1, 0, 0}, {0, 1, 0}, {0, 0, -1}},
{{1, 0, 0}, {0, 1, 0}, {0, 0, 1}}}));
// A right-angled triangle with side p.
>From bd48fe7859d80472350f97530949fdfcb1794b80 Mon Sep 17 00:00:00 2001
From: Abhinav271828 <abhinav.m at research.iiit.ac.in>
Date: Fri, 26 Jan 2024 22:37:00 +0530
Subject: [PATCH 36/41] Use parseRelationFromSet
---
.../Analysis/Presburger/BarvinokTest.cpp | 38 ++++++-------------
1 file changed, 11 insertions(+), 27 deletions(-)
diff --git a/mlir/unittests/Analysis/Presburger/BarvinokTest.cpp b/mlir/unittests/Analysis/Presburger/BarvinokTest.cpp
index 6799dd37caf1f..58669ebf71a48 100644
--- a/mlir/unittests/Analysis/Presburger/BarvinokTest.cpp
+++ b/mlir/unittests/Analysis/Presburger/BarvinokTest.cpp
@@ -1,5 +1,6 @@
#include "mlir/Analysis/Presburger/Barvinok.h"
#include "./Utils.h"
+#include "Parser.h"
#include <gmock/gmock.h>
#include <gtest/gtest.h>
@@ -243,16 +244,10 @@ TEST(BarvinokTest, computeNumTermsCone) {
/// hand.
TEST(BarvinokTest, computeNumTermsPolytope) {
// A cube of side 1.
- 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));
+ PolyhedronH poly =
+ parseRelationFromSet("(x, y, z) : (x >= 0, y >= 0, z >= 0, -x + 1 >= 0, "
+ "-y + 1 >= 0, -z + 1 >= 0)",
+ 0);
std::vector<std::pair<PresburgerSet, GeneratingFunction>> count =
computePolytopeGeneratingFunction(poly);
@@ -279,10 +274,8 @@ TEST(BarvinokTest, computeNumTermsPolytope) {
{{1, 0, 0}, {0, 1, 0}, {0, 0, 1}}}));
// A right-angled triangle with side p.
- ineqs = makeIntMatrix(3, 4, {{1, 0, 0, 0}, {0, 1, 0, 0}, {-1, -1, 1, 0}});
- poly = defineHRep(2, 1);
- for (unsigned i = 0; i < 3; i++)
- poly.addInequality(ineqs.getRow(i));
+ poly =
+ parseRelationFromSet("(x, y)[N] : (x >= 0, y >= 0, -x - y + N >= 0)", 0);
count = computePolytopeGeneratingFunction(poly);
// There is only one chamber: p ≥ 0
@@ -298,19 +291,10 @@ TEST(BarvinokTest, computeNumTermsPolytope) {
{{{-1, 1}, {-1, 0}}, {{1, -1}, {0, -1}}, {{1, 0}, {0, 1}}}));
// Cartesian product of a cube with side M and a right triangle with side N.
- ineqs = makeIntMatrix(9, 8,
- {{1, 0, 0, 0, 0, 1, 0, 0},
- {0, 1, 0, 0, 0, 1, 0, 0},
- {0, 0, 1, 0, 0, 1, 0, 0},
- {-1, 0, 0, 0, 0, 1, 0, 0},
- {0, -1, 0, 0, 0, 1, 0, 0},
- {0, 0, -1, 0, 0, 1, 0, 0},
- {0, 0, 0, 1, 0, 0, 0, 0},
- {0, 0, 0, 0, 1, 0, 0, 0},
- {0, 0, 0, -1, -1, 0, 1, 0}});
- poly = defineHRep(5, 2);
- for (unsigned i = 0; i < 9; i++)
- poly.addInequality(ineqs.getRow(i));
+ poly = parseRelationFromSet(
+ "(x, y, z, w, a)[M, N] : (x >= 0, y >= 0, z >= 0, -x + M >= 0, -y + M >= "
+ "0, -z + M >= 0, w >= 0, a >= 0, -w - a + N >= 0)",
+ 0);
count = computePolytopeGeneratingFunction(poly);
>From 36f035118508756a2005c3299f03e0a24ec74b07 Mon Sep 17 00:00:00 2001
From: Abhinav271828 <abhinav.m at research.iiit.ac.in>
Date: Mon, 29 Jan 2024 23:55:19 +0530
Subject: [PATCH 37/41] Move declaration
---
mlir/lib/Analysis/Presburger/Barvinok.cpp | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
diff --git a/mlir/lib/Analysis/Presburger/Barvinok.cpp b/mlir/lib/Analysis/Presburger/Barvinok.cpp
index fab8ac28af60e..558b3e72fe5ef 100644
--- a/mlir/lib/Analysis/Presburger/Barvinok.cpp
+++ b/mlir/lib/Analysis/Presburger/Barvinok.cpp
@@ -406,7 +406,6 @@ mlir::presburger::detail::computePolytopeGeneratingFunction(
// Now, we compute the generating function at this vertex.
// We collect the inequalities corresponding to each vertex to compute
// the tangent cone at that vertex.
- SmallVector<MPInt> ineq(numVars + 1);
// We only need the coefficients of the variables (NOT the parameters)
// as the generating function only depends on these.
@@ -414,6 +413,7 @@ mlir::presburger::detail::computePolytopeGeneratingFunction(
// constant terms zero.
ConeH tangentCone = defineHRep(numVars);
for (unsigned j = 0, e = subset.getNumRows(); j < e; ++j) {
+ SmallVector<MPInt> ineq(numVars + 1);
for (unsigned k = 0; k < numVars; ++k)
ineq[k] = subset(j, k);
tangentCone.addInequality(ineq);
>From 9d5d07b6e94a34a3864827e1b8d68dcf1b103c31 Mon Sep 17 00:00:00 2001
From: Abhinav271828 <abhinav.m at research.iiit.ac.in>
Date: Tue, 30 Jan 2024 00:03:42 +0530
Subject: [PATCH 38/41] Check for previously computed vertices
---
.../include/mlir/Analysis/Presburger/Matrix.h | 2 ++
mlir/lib/Analysis/Presburger/Barvinok.cpp | 21 ++++---------------
mlir/lib/Analysis/Presburger/Matrix.cpp | 18 ++++++++++++++++
3 files changed, 24 insertions(+), 17 deletions(-)
diff --git a/mlir/include/mlir/Analysis/Presburger/Matrix.h b/mlir/include/mlir/Analysis/Presburger/Matrix.h
index 15b39a4fd927f..2389be5fbf267 100644
--- a/mlir/include/mlir/Analysis/Presburger/Matrix.h
+++ b/mlir/include/mlir/Analysis/Presburger/Matrix.h
@@ -74,6 +74,8 @@ class Matrix {
T operator()(unsigned row, unsigned column) const { return at(row, column); }
+ bool operator==(const Matrix<T> &m) const;
+
/// Swap the given columns.
void swapColumns(unsigned column, unsigned otherColumn);
diff --git a/mlir/lib/Analysis/Presburger/Barvinok.cpp b/mlir/lib/Analysis/Presburger/Barvinok.cpp
index 558b3e72fe5ef..2fbf993f605e7 100644
--- a/mlir/lib/Analysis/Presburger/Barvinok.cpp
+++ b/mlir/lib/Analysis/Presburger/Barvinok.cpp
@@ -355,7 +355,10 @@ mlir::presburger::detail::computePolytopeGeneratingFunction(
if (!vertex)
continue;
- // If this subset corresponds to a vertex, store it.
+ if (std::find(vertices.begin(), vertices.end(), vertex) != vertices.end())
+ continue;
+ // If this subset corresponds to a vertex that has not been considered,
+ // store it.
vertices.push_back(*vertex);
// Let the current vertex be [X | y], where
@@ -381,22 +384,6 @@ mlir::presburger::detail::computePolytopeGeneratingFunction(
activeRegion.addToRow(i, b2c2.getRow(i), 1);
}
- // If any row of activeRegion is all-zero, that means that the
- // corresponding inequality in `remaining` is *also* satisfied by the
- // vertex for any values of the parameters. Thus we remove it from
- // activeRegion and add that inequality to `subset`.
- // Note that if the row is not all-zero, it is still possible for the
- // inequality to be satisfied by the corresponding vertex *in some subset of
- // the parameter space*. However, since this subset is defined by an
- // equality, it is not full-dimensional; here and in the chamber
- // decomposition, we ignore full-dimensional chambers.
- for (int i = activeRegion.getNumRows() - 1; i >= 0; --i) {
- if (!isRangeZero(activeRegion.getRow(i)))
- continue;
- subset.appendExtraRow(remainder.getRow(i));
- activeRegion.removeRow(i);
- }
-
// We convert the representation of the active region to an integers-only
// form so as to store it as a PresburgerSet.
IntegerPolyhedron activeRegionRel(
diff --git a/mlir/lib/Analysis/Presburger/Matrix.cpp b/mlir/lib/Analysis/Presburger/Matrix.cpp
index 3b251a1ca2266..0cbaacdeb2744 100644
--- a/mlir/lib/Analysis/Presburger/Matrix.cpp
+++ b/mlir/lib/Analysis/Presburger/Matrix.cpp
@@ -29,6 +29,24 @@ Matrix<T>::Matrix(unsigned rows, unsigned columns, unsigned reservedRows,
data.reserve(std::max(nRows, reservedRows) * nReservedColumns);
}
+template <typename T>
+bool Matrix<T>::operator==(const Matrix<T> &m) const {
+ if (nRows != m.getNumRows())
+ return false;
+ if (nColumns != m.getNumColumns())
+ return false;
+
+ for (unsigned i = 0; i < nRows; i++) {
+ for (unsigned j = 0; j < nColumns; j++) {
+ if (at(i, j) == m.at(i, j))
+ continue;
+ return false;
+ }
+ }
+
+ return true;
+}
+
template <typename T>
Matrix<T> Matrix<T>::identity(unsigned dimension) {
Matrix matrix(dimension, dimension);
>From 19a812a323e956644b70f992522b8e8059926797 Mon Sep 17 00:00:00 2001
From: Abhinav271828 <abhinav.m at research.iiit.ac.in>
Date: Tue, 30 Jan 2024 22:20:48 +0530
Subject: [PATCH 39/41] Fixes
---
.../mlir/Analysis/Presburger/Barvinok.h | 17 ++---
.../Analysis/Presburger/IntegerRelation.h | 8 +++
.../include/mlir/Analysis/Presburger/Matrix.h | 2 +-
.../Analysis/Presburger/PresburgerRelation.h | 6 +-
mlir/lib/Analysis/Presburger/Barvinok.cpp | 24 ++++---
mlir/lib/Analysis/Presburger/Matrix.cpp | 2 +-
.../Analysis/Presburger/BarvinokTest.cpp | 66 +++++++++----------
7 files changed, 63 insertions(+), 62 deletions(-)
diff --git a/mlir/include/mlir/Analysis/Presburger/Barvinok.h b/mlir/include/mlir/Analysis/Presburger/Barvinok.h
index a80aa9bc2f9b2..7c4c792eb32a1 100644
--- a/mlir/include/mlir/Analysis/Presburger/Barvinok.h
+++ b/mlir/include/mlir/Analysis/Presburger/Barvinok.h
@@ -95,16 +95,13 @@ GeneratingFunction computeUnimodularConeGeneratingFunction(ParamPoint vertex,
/// If there is no solution, return null.
std::optional<ParamPoint> solveParametricEquations(FracMatrix equations);
-/// Given a list of possibly intersecting regions (PresburgerSet) 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]}.
-/// The ith relation corresponds to the activity region of the ith vertex set.
-/// Note that here, by disjoint, we mean that the intersection is not
-/// full-dimensional.
+/// Given a list of possibly intersecting regions (PresburgerSet) and the
+/// generating functions active in each region, produce a pairwise disjoint
+/// list of regions (chambers) and identify the generating function of the
+/// polytope in each chamber.
+/// The returned list partitions the universe into parts depending on which
+/// subset of GFs is active there, and gives the sum of active GFs for each
+/// part. Lower-dimensional partitions are ignored.
std::vector<std::pair<PresburgerSet, GeneratingFunction>>
computeChamberDecomposition(
unsigned numSymbols, ArrayRef<std::pair<PresburgerSet, GeneratingFunction>>
diff --git a/mlir/include/mlir/Analysis/Presburger/IntegerRelation.h b/mlir/include/mlir/Analysis/Presburger/IntegerRelation.h
index b4159b490381e..0fc786ee88a32 100644
--- a/mlir/include/mlir/Analysis/Presburger/IntegerRelation.h
+++ b/mlir/include/mlir/Analysis/Presburger/IntegerRelation.h
@@ -891,6 +891,14 @@ class IntegerPolyhedron : public IntegerRelation {
addInequality(inequalities.getRow(i));
}
+ explicit IntegerPolyhedron(const PresburgerSpace &space,
+ FracMatrix inequalities)
+ : IntegerPolyhedron(space) {
+ IntMatrix ineqsNormalized = inequalities.normalizeRows();
+ for (unsigned i = 0, e = inequalities.getNumRows(); i < e; i++)
+ addInequality(ineqsNormalized.getRow(i));
+ }
+
/// Construct a set from an IntegerRelation. The relation should have
/// no domain vars.
explicit IntegerPolyhedron(const IntegerRelation &rel)
diff --git a/mlir/include/mlir/Analysis/Presburger/Matrix.h b/mlir/include/mlir/Analysis/Presburger/Matrix.h
index 2389be5fbf267..7f2e7c2e5b84d 100644
--- a/mlir/include/mlir/Analysis/Presburger/Matrix.h
+++ b/mlir/include/mlir/Analysis/Presburger/Matrix.h
@@ -203,7 +203,7 @@ class Matrix {
/// 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::vector<int> indicator);
+ std::pair<Matrix<T>, Matrix<T>> splitByBitset(ArrayRef<int> indicator);
/// Print the matrix.
void print(raw_ostream &os) const;
diff --git a/mlir/include/mlir/Analysis/Presburger/PresburgerRelation.h b/mlir/include/mlir/Analysis/Presburger/PresburgerRelation.h
index 7d28176f477a3..08b5007cc67ea 100644
--- a/mlir/include/mlir/Analysis/Presburger/PresburgerRelation.h
+++ b/mlir/include/mlir/Analysis/Presburger/PresburgerRelation.h
@@ -217,9 +217,9 @@ class PresburgerRelation {
/// redundencies.
PresburgerRelation simplify() const;
- /// Return true if any of the disjuncts is full-dimensional.
- /// We find if a disjunct is full-dimensional by checking if it is flat
- /// along the dimension of any of its inequalities.
+ /// We consider a disjunct to be full-dimensional if it is not flat along the
+ /// dimension of any of its inequalities.
+ /// A PresburgerRelation is full-dimensional if any of its disjuncts is.
bool isFullDim() const;
/// Print the set's internal state.
diff --git a/mlir/lib/Analysis/Presburger/Barvinok.cpp b/mlir/lib/Analysis/Presburger/Barvinok.cpp
index 2fbf993f605e7..fafbc0a12d294 100644
--- a/mlir/lib/Analysis/Presburger/Barvinok.cpp
+++ b/mlir/lib/Analysis/Presburger/Barvinok.cpp
@@ -224,16 +224,15 @@ mlir::presburger::detail::solveParametricEquations(FracMatrix equations) {
/// This is an implementation of the Clauss-Loechner algorithm for chamber
/// decomposition.
-/// We maintain a list of pairwise disjoint chambers and their generating
-/// functions. We iterate over the list of regions, each time adding the
-/// generating function to the chambers where the corresponding vertex is
-/// active and separating the chambers where it is not.
+/// We maintain a list of pairwise disjoint chambers and the generating
+/// functions corresponding to each one. We iterate over the list of regions,
+/// each time adding the current region's generating function to the chambers
+/// where it is active and separating the chambers where it is not.
///
-/// Given the region each vertex is active in, for each subset of vertices,
-/// the region that precisely this subset is in, is the intersection of the
-/// regions that these are active in, intersected with the complements of the
-/// remaining regions. The generating function for this region is the sum of
-/// generating functions for the vertices' tangent cones.
+/// Given the region each generating function is active in, for each subset of
+/// generating functions the region that (the sum of) precisely this subset is
+/// in, is the intersection of the regions that these are active in,
+/// intersected with the complements of the remaining regions.
std::vector<std::pair<PresburgerSet, GeneratingFunction>>
mlir::presburger::detail::computeChamberDecomposition(
unsigned numSymbols, ArrayRef<std::pair<PresburgerSet, GeneratingFunction>>
@@ -329,7 +328,7 @@ mlir::presburger::detail::computePolytopeGeneratingFunction(
// For a given permutation, we consider a subset which contains
// the i'th inequality if the i'th bit in the bitset is 1.
// We start with the permutation that takes the last numVars inequalities.
- std::vector<int> indicator(numIneqs);
+ SmallVector<int> indicator(numIneqs);
for (unsigned i = numIneqs - numVars; i < numIneqs; ++i)
indicator[i] = 1;
@@ -368,7 +367,7 @@ mlir::presburger::detail::computePolytopeGeneratingFunction(
// 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].
+ // remaining = [A2 | B2 | c2].
// Thus, the coefficients of the parameters after substitution become
// (A2 • X + B2)
// and the constant terms become
@@ -387,8 +386,7 @@ mlir::presburger::detail::computePolytopeGeneratingFunction(
// We convert the representation of the active region to an integers-only
// form so as to store it as a PresburgerSet.
IntegerPolyhedron activeRegionRel(
- PresburgerSpace::getRelationSpace(0, numSymbols, 0, 0),
- activeRegion.normalizeRows());
+ PresburgerSpace::getRelationSpace(0, numSymbols, 0, 0), activeRegion);
// Now, we compute the generating function at this vertex.
// We collect the inequalities corresponding to each vertex to compute
diff --git a/mlir/lib/Analysis/Presburger/Matrix.cpp b/mlir/lib/Analysis/Presburger/Matrix.cpp
index 0cbaacdeb2744..484fb05adab05 100644
--- a/mlir/lib/Analysis/Presburger/Matrix.cpp
+++ b/mlir/lib/Analysis/Presburger/Matrix.cpp
@@ -413,7 +413,7 @@ void Matrix<T>::print(raw_ostream &os) const {
/// 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::vector<int> indicator) {
+Matrix<T>::splitByBitset(ArrayRef<int> indicator) {
Matrix<T> rowsForOne(0, nColumns), rowsForZero(0, nColumns);
for (unsigned i = 0; i < nRows; i++) {
if (indicator[i] == 1)
diff --git a/mlir/unittests/Analysis/Presburger/BarvinokTest.cpp b/mlir/unittests/Analysis/Presburger/BarvinokTest.cpp
index 58669ebf71a48..fedfd5a55c9a9 100644
--- a/mlir/unittests/Analysis/Presburger/BarvinokTest.cpp
+++ b/mlir/unittests/Analysis/Presburger/BarvinokTest.cpp
@@ -52,40 +52,38 @@ TEST(BarvinokTest, getIndex) {
// (s.t. the cones are unimodular) and the generating functions
// are computed. We check that the results contain the correct
// matrices.
-// TEST(BarvinokTest, unimodularConeGeneratingFunction) {
-// ConeH cone = defineHRep(2);
-// cone.addInequality({0, -1, 0});
-// cone.addInequality({-1, -2, 0});
-//
-// ParamPoint vertex =
-// makeFracMatrix(2, 3, {{2, 2, 0}, {-1, -Fraction(1, 2), 1}});
-//
-// GeneratingFunction gf =
-// computeUnimodularConeGeneratingFunction(vertex, 1, cone);
-//
-// EXPECT_EQ_REPR_GENERATINGFUNCTION(
-// gf, GeneratingFunction(
-// 2, {1},
-// {makeFracMatrix(3, 2, {{-1, 0}, {-Fraction(1, 2), 1}, {1,
-// 2}})},
-// {{{2, -1}, {-1, 0}}}));
-//
-// cone = defineHRep(3);
-// cone.addInequality({7, 1, 6, 0});
-// cone.addInequality({9, 1, 7, 0});
-// cone.addInequality({8, -1, 1, 0});
-//
-// vertex = makeFracMatrix(3, 2, {{5, 2}, {6, 2}, {7, 1}});
-//
-// gf = computeUnimodularConeGeneratingFunction(vertex, 1, cone);
-//
-// EXPECT_EQ_REPR_GENERATINGFUNCTION(
-// gf,
-// GeneratingFunction(
-// 1, {1}, {makeFracMatrix(2, 3, {{-83, -100, -41}, {-22, -27,
-// -15}})},
-// {{{8, 47, -17}, {-7, -41, 15}, {1, 5, -2}}}));
-// }
+TEST(BarvinokTest, unimodularConeGeneratingFunction) {
+ ConeH cone = defineHRep(2);
+ cone.addInequality({0, -1, 0});
+ cone.addInequality({-1, -2, 0});
+
+ ParamPoint vertex =
+ makeFracMatrix(2, 3, {{2, 2, 0}, {-1, -Fraction(1, 2), 1}});
+
+ GeneratingFunction gf =
+ computeUnimodularConeGeneratingFunction(vertex, 1, cone);
+
+ EXPECT_EQ_REPR_GENERATINGFUNCTION(
+ gf, GeneratingFunction(
+ 2, {1},
+ {makeFracMatrix(3, 2, {{-1, 0}, {-Fraction(1, 2), 1}, {1, 2}})},
+ {{{2, -1}, {-1, 0}}}));
+
+ cone = defineHRep(3);
+ cone.addInequality({7, 1, 6, 0});
+ cone.addInequality({9, 1, 7, 0});
+ cone.addInequality({8, -1, 1, 0});
+
+ vertex = makeFracMatrix(3, 2, {{5, 2}, {6, 2}, {7, 1}});
+
+ gf = computeUnimodularConeGeneratingFunction(vertex, 1, cone);
+
+ EXPECT_EQ_REPR_GENERATINGFUNCTION(
+ gf,
+ GeneratingFunction(
+ 1, {1}, {makeFracMatrix(2, 3, {{-83, -100, -41}, {-22, -27, -15}})},
+ {{{8, 47, -17}, {-7, -41, 15}, {1, 5, -2}}}));
+}
// The following vectors are randomly generated.
// We then check that the output of the function has non-zero
>From 709aba0ee331ec880b39a8fb20da0c842eb88dfa Mon Sep 17 00:00:00 2001
From: Abhinav271828 <abhinav.m at research.iiit.ac.in>
Date: Tue, 30 Jan 2024 22:21:59 +0530
Subject: [PATCH 40/41] doc
---
mlir/lib/Analysis/Presburger/Barvinok.cpp | 4 ++--
1 file changed, 2 insertions(+), 2 deletions(-)
diff --git a/mlir/lib/Analysis/Presburger/Barvinok.cpp b/mlir/lib/Analysis/Presburger/Barvinok.cpp
index fafbc0a12d294..9f57cfed2d583 100644
--- a/mlir/lib/Analysis/Presburger/Barvinok.cpp
+++ b/mlir/lib/Analysis/Presburger/Barvinok.cpp
@@ -366,8 +366,8 @@ mlir::presburger::detail::computePolytopeGeneratingFunction(
// 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.,
- // remaining = [A2 | B2 | c2].
+ // polytope (those which were not collected into `subset`), i.e., into the
+ // matrix [A2 | B2 | c2].
// Thus, the coefficients of the parameters after substitution become
// (A2 • X + B2)
// and the constant terms become
>From c8de2377ce7c7a22de2bff438eafc1ca094a42db Mon Sep 17 00:00:00 2001
From: Abhinav271828 <abhinav.m at research.iiit.ac.in>
Date: Tue, 30 Jan 2024 23:18:53 +0530
Subject: [PATCH 41/41] Fix
---
.../mlir/Analysis/Presburger/Barvinok.h | 1 +
.../Analysis/Presburger/PresburgerRelation.h | 5 ++---
mlir/lib/Analysis/Presburger/Barvinok.cpp | 21 +++++++++++++++----
mlir/lib/Analysis/Presburger/Matrix.cpp | 10 ++++-----
4 files changed, 24 insertions(+), 13 deletions(-)
diff --git a/mlir/include/mlir/Analysis/Presburger/Barvinok.h b/mlir/include/mlir/Analysis/Presburger/Barvinok.h
index 7c4c792eb32a1..65e116e0aa1fb 100644
--- a/mlir/include/mlir/Analysis/Presburger/Barvinok.h
+++ b/mlir/include/mlir/Analysis/Presburger/Barvinok.h
@@ -99,6 +99,7 @@ std::optional<ParamPoint> solveParametricEquations(FracMatrix equations);
/// generating functions active in each region, produce a pairwise disjoint
/// list of regions (chambers) and identify the generating function of the
/// polytope in each chamber.
+///
/// The returned list partitions the universe into parts depending on which
/// subset of GFs is active there, and gives the sum of active GFs for each
/// part. Lower-dimensional partitions are ignored.
diff --git a/mlir/include/mlir/Analysis/Presburger/PresburgerRelation.h b/mlir/include/mlir/Analysis/Presburger/PresburgerRelation.h
index 08b5007cc67ea..9634df6d58a1a 100644
--- a/mlir/include/mlir/Analysis/Presburger/PresburgerRelation.h
+++ b/mlir/include/mlir/Analysis/Presburger/PresburgerRelation.h
@@ -217,9 +217,8 @@ class PresburgerRelation {
/// redundencies.
PresburgerRelation simplify() const;
- /// We consider a disjunct to be full-dimensional if it is not flat along the
- /// dimension of any of its inequalities.
- /// A PresburgerRelation is full-dimensional if any of its disjuncts is.
+ /// Return whether the given PresburgerRelation is full-dimensional. By full-
+ /// dimensional we mean that it is not flat along any dimension.
bool isFullDim() const;
/// Print the set's internal state.
diff --git a/mlir/lib/Analysis/Presburger/Barvinok.cpp b/mlir/lib/Analysis/Presburger/Barvinok.cpp
index 9f57cfed2d583..42d4bccf47ca6 100644
--- a/mlir/lib/Analysis/Presburger/Barvinok.cpp
+++ b/mlir/lib/Analysis/Presburger/Barvinok.cpp
@@ -155,6 +155,7 @@ mlir::presburger::detail::computeUnimodularConeGeneratingFunction(
/// 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.
@@ -211,9 +212,11 @@ mlir::presburger::detail::solveParametricEquations(FracMatrix equations) {
// 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 =
equations.getSubMatrix(/*fromRow=*/0, /*toRow=*/d - 1,
@@ -224,6 +227,7 @@ mlir::presburger::detail::solveParametricEquations(FracMatrix equations) {
/// This is an implementation of the Clauss-Loechner algorithm for chamber
/// decomposition.
+///
/// We maintain a list of pairwise disjoint chambers and the generating
/// functions corresponding to each one. We iterate over the list of regions,
/// each time adding the current region's generating function to the chambers
@@ -246,18 +250,21 @@ mlir::presburger::detail::computeChamberDecomposition(
GeneratingFunction(numSymbols, {}, {}, {})}};
// We iterate over the region list.
+ //
// For each activity region R_j (corresponding to a vertex v_j, whose
// generating function is gf_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 (if it is full-dimensional), where the
// generating function is gf_i + gf_j.
// 2. the difference R_i - R_j (if it is full-dimensional), where v_j is
// inactive and the generating function is gf_i.
-
+ //
// At each step, we define a new chamber list after considering vertex v_j,
// and its generating function, replacing and appending chambers as
// discussed above.
+ //
// The loop has the invariant that the union over all the chambers gives the
// universe at every step (modulo lower-dimensional spaces).
for (const auto &[region, generatingFunction] :
@@ -367,14 +374,18 @@ mlir::presburger::detail::computePolytopeGeneratingFunction(
// 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., into the
- // matrix [A2 | B2 | c2].
+ // inequalities [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]
+ //
+ // 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].
FracMatrix activeRegion(numIneqs - numVars, numSymbols + 1);
@@ -405,6 +416,7 @@ mlir::presburger::detail::computePolytopeGeneratingFunction(
}
// We assume that the tangent cone is unimodular, so there is no need
// to decompose it.
+ //
// In the general case, the unimodular decomposition may have several
// cones.
GeneratingFunction vertexGf(numSymbols, {}, {}, {});
@@ -415,7 +427,7 @@ mlir::presburger::detail::computePolytopeGeneratingFunction(
computeUnimodularConeGeneratingFunction(*vertex, sign, cone);
}
// We store the vertex we computed with the generating function of its
- // tangent cones.
+ // tangent cone.
regionsAndGeneratingFunctions.emplace_back(PresburgerSet(activeRegionRel),
vertexGf);
} while (std::next_permutation(indicator.begin(), indicator.end()));
@@ -425,6 +437,7 @@ mlir::presburger::detail::computePolytopeGeneratingFunction(
// property that no two of them have a full-dimensional intersection, i.e.,
// they may share "facets" or "edges", but their intersection can only have
// up to numVars - 1 dimensions.
+ //
// In each chamber, we sum up the generating functions of the active vertices
// to find the generating function of the polytope.
return computeChamberDecomposition(numSymbols, regionsAndGeneratingFunctions);
diff --git a/mlir/lib/Analysis/Presburger/Matrix.cpp b/mlir/lib/Analysis/Presburger/Matrix.cpp
index 484fb05adab05..4cb6e6b16bc87 100644
--- a/mlir/lib/Analysis/Presburger/Matrix.cpp
+++ b/mlir/lib/Analysis/Presburger/Matrix.cpp
@@ -29,6 +29,8 @@ Matrix<T>::Matrix(unsigned rows, unsigned columns, unsigned reservedRows,
data.reserve(std::max(nRows, reservedRows) * nReservedColumns);
}
+/// We cannot use the default implementation of operator== as it compares
+/// fields like `reservedColumns` etc., which are not part of the data.
template <typename T>
bool Matrix<T>::operator==(const Matrix<T> &m) const {
if (nRows != m.getNumRows())
@@ -36,13 +38,9 @@ bool Matrix<T>::operator==(const Matrix<T> &m) const {
if (nColumns != m.getNumColumns())
return false;
- for (unsigned i = 0; i < nRows; i++) {
- for (unsigned j = 0; j < nColumns; j++) {
- if (at(i, j) == m.at(i, j))
- continue;
+ for (unsigned i = 0; i < nRows; i++)
+ if (getRow(i) != m.getRow(i))
return false;
- }
- }
return true;
}
More information about the Mlir-commits
mailing list