[clang-tools-extra] [mlir] [llvm] [compiler-rt] [flang] [clang] [MLIR][Presburger] Function to compute the generating function corresponding to a unimodular cone (PR #77199)

via cfe-commits cfe-commits at lists.llvm.org
Sat Jan 6 06:06:56 PST 2024


llvmbot wrote:


<!--LLVM PR SUMMARY COMMENT-->

@llvm/pr-subscribers-mlir-presburger

Author: None (Abhinav271828)

<details>
<summary>Changes</summary>

We implement a function to compute the generating function for a signed unimodular cone with a parametric vertex.
These functions, summed over tangential cones, provide the generating function corresponding to a parametric polyhedron.

---
Full diff: https://github.com/llvm/llvm-project/pull/77199.diff


8 Files Affected:

- (added) mlir/include/mlir/Analysis/Presburger/Barvinok.h (+90) 
- (renamed) mlir/include/mlir/Analysis/Presburger/GeneratingFunction.h (+4-1) 
- (added) mlir/lib/Analysis/Presburger/Barvinok.cpp (+134) 
- (modified) mlir/lib/Analysis/Presburger/CMakeLists.txt (+1) 
- (added) mlir/unittests/Analysis/Presburger/BarvinokTest.cpp (+48) 
- (modified) mlir/unittests/Analysis/Presburger/CMakeLists.txt (+2) 
- (added) mlir/unittests/Analysis/Presburger/GeneratingFunctionTest.cpp (+39) 
- (modified) mlir/unittests/Analysis/Presburger/Utils.h (+35-1) 


``````````diff
diff --git a/mlir/include/mlir/Analysis/Presburger/Barvinok.h b/mlir/include/mlir/Analysis/Presburger/Barvinok.h
new file mode 100644
index 00000000000000..93b29e2d718e59
--- /dev/null
+++ b/mlir/include/mlir/Analysis/Presburger/Barvinok.h
@@ -0,0 +1,90 @@
+//===- Barvinok.h - Barvinok's Algorithm ------------------------*- C++ -*-===//
+//
+// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
+// See https://llvm.org/LICENSE.txt for license information.
+// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
+//
+//===----------------------------------------------------------------------===//
+//
+// Implementation of Barvinok's algorithm and related utility functions.
+// Currently a work in progress.
+// These include functions to manipulate cones (define a cone object, get its
+// dual, and find its index).
+//
+// The implementation is based on:
+// 1. Barvinok, Alexander, and James E. Pommersheim. "An algorithmic theory of
+//    lattice points in polyhedra." New perspectives in algebraic combinatorics
+//    38 (1999): 91-147.
+// 2. Verdoolaege, Sven, et al. "Counting integer points in parametric
+//    polytopes using Barvinok's rational functions." Algorithmica 48 (2007):
+//    37-66.
+//
+//===----------------------------------------------------------------------===//
+
+#ifndef MLIR_ANALYSIS_PRESBURGER_BARVINOK_H
+#define MLIR_ANALYSIS_PRESBURGER_BARVINOK_H
+
+#include "mlir/Analysis/Presburger/GeneratingFunction.h"
+#include "mlir/Analysis/Presburger/IntegerRelation.h"
+#include "mlir/Analysis/Presburger/Matrix.h"
+#include <optional>
+
+namespace mlir {
+namespace presburger {
+namespace detail {
+
+/// A polyhedron in H-representation is a set of inequalities
+/// in d variables with integer coefficients.
+using PolyhedronH = IntegerRelation;
+
+/// A polyhedron in V-representation is a set of rays and points, i.e.,
+/// vectors, stored as rows of a matrix.
+using PolyhedronV = IntMatrix;
+
+/// A cone in either representation is a special case of
+/// a polyhedron in that representation.
+using ConeH = PolyhedronH;
+using ConeV = PolyhedronV;
+
+inline ConeH defineHRep(int numVars) {
+  // 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.
+  // Once the cone is defined, we use `addInequality()` to set inequalities.
+  return ConeH(PresburgerSpace::getSetSpace(/*numDims=*/numVars,
+                                            /*numSymbols=*/0,
+                                            /*numLocals=*/0));
+}
+
+/// Get the index of a cone, i.e., the volume of the parallelepiped
+/// spanned by its generators, which is equal to the number of integer
+/// points in its fundamental parallelepiped.
+/// If the index is 1, the cone is unimodular.
+/// Barvinok, A., and J. E. Pommersheim. "An algorithmic theory of lattice
+/// points in polyhedra." p. 107 If it has more rays than the dimension, return
+/// 0.
+MPInt getIndex(ConeV cone);
+
+/// Given a cone in H-representation, return its dual. The dual cone is in
+/// V-representation.
+/// This assumes that the input is pointed at the origin; it assert-fails
+/// otherwise.
+ConeV getDual(ConeH cone);
+
+/// Given a cone in V-representation, return its dual. The dual cone is in
+/// H-representation.
+/// The returned cone is pointed at the origin.
+ConeH getDual(ConeV cone);
+
+/// Compute the generating function for a unimodular cone.
+/// It assert-fails if the input cone is not unimodular.
+GeneratingFunction unimodularConeGeneratingFunction(ParamPoint vertex, int sign,
+                                                    ConeH cone);
+
+} // namespace detail
+} // namespace presburger
+} // namespace mlir
+
+#endif // MLIR_ANALYSIS_PRESBURGER_BARVINOK_H
diff --git a/mlir/lib/Analysis/Presburger/GeneratingFunction.h b/mlir/include/mlir/Analysis/Presburger/GeneratingFunction.h
similarity index 96%
rename from mlir/lib/Analysis/Presburger/GeneratingFunction.h
rename to mlir/include/mlir/Analysis/Presburger/GeneratingFunction.h
index f7deba921ea51e..4dd692c251563b 100644
--- a/mlir/lib/Analysis/Presburger/GeneratingFunction.h
+++ b/mlir/include/mlir/Analysis/Presburger/GeneratingFunction.h
@@ -19,6 +19,7 @@
 
 namespace mlir {
 namespace presburger {
+namespace detail {
 
 // A parametric point is a vector, each of whose elements
 // is an affine function of n parameters. Each row
@@ -83,7 +84,8 @@ class GeneratingFunction {
     std::vector<std::vector<Point>> sumDenominators = denominators;
     sumDenominators.insert(sumDenominators.end(), gf.denominators.begin(),
                            gf.denominators.end());
-    return GeneratingFunction(0, sumSigns, sumNumerators, sumDenominators);
+    return GeneratingFunction(numParam, sumSigns, sumNumerators,
+                              sumDenominators);
   }
 
   llvm::raw_ostream &print(llvm::raw_ostream &os) const {
@@ -128,6 +130,7 @@ class GeneratingFunction {
   std::vector<std::vector<Point>> denominators;
 };
 
+} // namespace detail
 } // namespace presburger
 } // namespace mlir
 
diff --git a/mlir/lib/Analysis/Presburger/Barvinok.cpp b/mlir/lib/Analysis/Presburger/Barvinok.cpp
new file mode 100644
index 00000000000000..31fafbda1129bf
--- /dev/null
+++ b/mlir/lib/Analysis/Presburger/Barvinok.cpp
@@ -0,0 +1,134 @@
+//===- Barvinok.cpp - Barvinok's Algorithm ----------------------*- C++ -*-===//
+//
+// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
+// See https://llvm.org/LICENSE.txt for license information.
+// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
+//
+//===----------------------------------------------------------------------===//
+
+#include "mlir/Analysis/Presburger/Barvinok.h"
+#include "mlir/Analysis/Presburger/GeneratingFunction.h"
+
+using namespace mlir;
+using namespace presburger;
+using namespace mlir::presburger::detail;
+
+/// Assuming that the input cone is pointed at the origin,
+/// converts it to its dual in V-representation.
+/// Essentially we just remove the all-zeroes constant column.
+ConeV mlir::presburger::detail::getDual(ConeH cone) {
+  unsigned numIneq = cone.getNumInequalities();
+  unsigned numVar = cone.getNumCols() - 1;
+  ConeV dual(numIneq, numVar, 0, 0);
+  // Assuming that an inequality of the form
+  // a1*x1 + ... + an*xn + b ≥ 0
+  // is represented as a row [a1, ..., an, b]
+  // and that b = 0.
+
+  for (unsigned i = 0; i < numIneq; ++i) {
+    assert(cone.atIneq(i, numVar) == 0 &&
+           "H-representation of cone is not centred at the origin!");
+    for (unsigned j = 0; j < numVar; ++j) {
+      dual.at(i, j) = cone.atIneq(i, j);
+    }
+  }
+
+  // Now dual is of the form [ [a1, ..., an] , ... ]
+  // which is the V-representation of the dual.
+  return dual;
+}
+
+/// Converts a cone in V-representation to the H-representation
+/// of its dual, pointed at the origin (not at the original vertex).
+/// Essentially adds a column consisting only of zeroes to the end.
+ConeH mlir::presburger::detail::getDual(ConeV cone) {
+  unsigned rows = cone.getNumRows();
+  unsigned columns = cone.getNumColumns();
+  ConeH dual = defineHRep(columns);
+  // Add a new column (for constants) at the end.
+  // This will be initialized to zero.
+  cone.insertColumn(columns);
+
+  for (unsigned i = 0; i < rows; ++i)
+    dual.addInequality(cone.getRow(i));
+
+  // Now dual is of the form [ [a1, ..., an, 0] , ... ]
+  // which is the H-representation of the dual.
+  return dual;
+}
+
+/// Find the index of a cone in V-representation.
+MPInt mlir::presburger::detail::getIndex(ConeV cone) {
+  if (cone.getNumRows() > cone.getNumColumns())
+    return MPInt(0);
+
+  return cone.determinant();
+}
+
+/// Compute the generating function for a unimodular cone.
+GeneratingFunction mlir::presburger::detail::unimodularConeGeneratingFunction(
+    ParamPoint vertex, int sign, ConeH cone) {
+  // `cone` is assumed to be unimodular.
+  assert(getIndex(getDual(cone)) == 1 && "input cone is not unimodular!");
+
+  unsigned numVar = cone.getNumVars();
+  unsigned numIneq = cone.getNumInequalities();
+
+  // Thus its ray matrix, U, is the inverse of the
+  // transpose of its inequality matrix, `cone`.
+  FracMatrix transp(numVar, numIneq);
+  for (unsigned i = 0; i < numVar; i++)
+    for (unsigned j = 0; j < numIneq; j++)
+      transp(j, i) = Fraction(cone.atIneq(i, j), 1);
+
+  FracMatrix generators(numVar, numIneq);
+  transp.determinant(&generators); // This is the U-matrix.
+
+  // The denominators of the generating function
+  // are given by the generators of the cone, i.e.,
+  // the rows of the matrix U.
+  std::vector<Point> denominator(numIneq);
+  ArrayRef<Fraction> row;
+  for (unsigned i = 0; i < numVar; i++) {
+    row = generators.getRow(i);
+    denominator[i] = Point(row);
+  }
+
+  // The vertex is v : [d, n+1].
+  // We need to find affine functions of parameters λi(p)
+  // such that v = Σ λi(p)*ui.
+  // The λi are given by the columns of Λ = v^T @ U^{-1} = v^T @ transp.
+  // Then the numerator will be Σ -floor(-λi(p))*u_i.
+  // Thus we store the numerator as the affine function -Λ,
+  // since the generators are already stored in the denominator.
+  // Note that the outer -1 will have to be accounted for, as it is not stored.
+  // See end for an example.
+
+  unsigned numColumns = vertex.getNumColumns();
+  unsigned numRows = vertex.getNumRows();
+  ParamPoint numerator(numColumns, numRows);
+  SmallVector<Fraction> ithCol(numRows);
+  for (unsigned i = 0; i < numColumns; i++) {
+    for (unsigned j = 0; j < numRows; j++)
+      ithCol[j] = vertex(j, i);
+    numerator.setRow(i, transp.preMultiplyWithRow(ithCol));
+    numerator.negateRow(i);
+  }
+
+  return GeneratingFunction(numColumns, SmallVector<int>(1, sign),
+                            std::vector({numerator}),
+                            std::vector({denominator}));
+
+  // Suppose the vertex is given by the matrix [ 2  2   0], with 2 params
+  //                                           [-1 -1/2 1]
+  // and the cone has H-representation [0  -1] => U-matrix [ 2 -1]
+  //                                   [-1 -2]             [-1  0]
+  // Therefore Λ will be given by [ 1    0 ] and the negation of this will be
+  // stored as the numerator.
+  //                              [ 1/2 -1 ]
+  //                              [ -1  -2 ]
+
+  // Algebraically, the numerator exponent is
+  // [ -2 ⌊ -N - M/2 +1 ⌋ + 1 ⌊ 0 +M +2 ⌋ ] -> first  COLUMN of U is [2, -1]
+  // [  1 ⌊ -N - M/2 +1 ⌋ + 0 ⌊ 0 +M +2 ⌋ ] -> second COLUMN of U is [-1, 0]
+}
diff --git a/mlir/lib/Analysis/Presburger/CMakeLists.txt b/mlir/lib/Analysis/Presburger/CMakeLists.txt
index e77e1623dae175..83d0514c9e7d17 100644
--- a/mlir/lib/Analysis/Presburger/CMakeLists.txt
+++ b/mlir/lib/Analysis/Presburger/CMakeLists.txt
@@ -1,4 +1,5 @@
 add_mlir_library(MLIRPresburger
+  Barvinok.cpp
   IntegerRelation.cpp
   LinearTransform.cpp
   Matrix.cpp
diff --git a/mlir/unittests/Analysis/Presburger/BarvinokTest.cpp b/mlir/unittests/Analysis/Presburger/BarvinokTest.cpp
new file mode 100644
index 00000000000000..b88baa6c6b48a4
--- /dev/null
+++ b/mlir/unittests/Analysis/Presburger/BarvinokTest.cpp
@@ -0,0 +1,48 @@
+#include "mlir/Analysis/Presburger/Barvinok.h"
+#include "./Utils.h"
+#include <gmock/gmock.h>
+#include <gtest/gtest.h>
+
+using namespace mlir;
+using namespace presburger;
+using namespace mlir::presburger::detail;
+
+/// The following are 3 randomly generated vectors with 4
+/// entries each and define a cone's H-representation
+/// using these numbers. We check that the dual contains
+/// the same numbers.
+/// We do the same in the reverse case.
+TEST(BarvinokTest, getDual) {
+  ConeH cone1 = defineHRep(4);
+  cone1.addInequality({1, 2, 3, 4, 0});
+  cone1.addInequality({3, 4, 2, 5, 0});
+  cone1.addInequality({6, 2, 6, 1, 0});
+
+  ConeV dual1 = getDual(cone1);
+
+  EXPECT_EQ_INT_MATRIX(
+      dual1, makeIntMatrix(3, 4, {{1, 2, 3, 4}, {3, 4, 2, 5}, {6, 2, 6, 1}}));
+
+  ConeV cone2 = makeIntMatrix(3, 4, {{3, 6, 1, 5}, {3, 1, 7, 2}, {9, 3, 2, 7}});
+
+  ConeH dual2 = getDual(cone2);
+
+  ConeH expected = defineHRep(4);
+  expected.addInequality({3, 6, 1, 5, 0});
+  expected.addInequality({3, 1, 7, 2, 0});
+  expected.addInequality({9, 3, 2, 7, 0});
+
+  EXPECT_TRUE(dual2.isEqual(expected));
+}
+
+/// We randomly generate a nxn matrix to use as a cone
+/// with n inequalities in n variables and check for
+/// the determinant being equal to the index.
+TEST(BarvinokTest, getIndex) {
+  ConeV cone = makeIntMatrix(3, 3, {{4, 2, 1}, {5, 2, 7}, {4, 1, 6}});
+  EXPECT_EQ(getIndex(cone), cone.determinant());
+
+  cone = makeIntMatrix(
+      4, 4, {{4, 2, 5, 1}, {4, 1, 3, 6}, {8, 2, 5, 6}, {5, 2, 5, 7}});
+  EXPECT_EQ(getIndex(cone), cone.determinant());
+}
diff --git a/mlir/unittests/Analysis/Presburger/CMakeLists.txt b/mlir/unittests/Analysis/Presburger/CMakeLists.txt
index e37133354e53ca..c98668f63fa5dc 100644
--- a/mlir/unittests/Analysis/Presburger/CMakeLists.txt
+++ b/mlir/unittests/Analysis/Presburger/CMakeLists.txt
@@ -1,5 +1,7 @@
 add_mlir_unittest(MLIRPresburgerTests
+  BarvinokTest.cpp
   FractionTest.cpp
+  GeneratingFunctionTest.cpp
   IntegerPolyhedronTest.cpp
   IntegerRelationTest.cpp
   LinearTransformTest.cpp
diff --git a/mlir/unittests/Analysis/Presburger/GeneratingFunctionTest.cpp b/mlir/unittests/Analysis/Presburger/GeneratingFunctionTest.cpp
new file mode 100644
index 00000000000000..5df1a5a0f64c05
--- /dev/null
+++ b/mlir/unittests/Analysis/Presburger/GeneratingFunctionTest.cpp
@@ -0,0 +1,39 @@
+//===- MatrixTest.cpp - Tests for QuasiPolynomial -------------------------===//
+//
+// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
+// See https://llvm.org/LICENSE.txt for license information.
+// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
+//
+//===----------------------------------------------------------------------===//
+
+#include "mlir/Analysis/Presburger/GeneratingFunction.h"
+#include "./Utils.h"
+#include <gmock/gmock.h>
+#include <gtest/gtest.h>
+
+using namespace mlir;
+using namespace presburger;
+using namespace mlir::presburger::detail;
+
+TEST(GeneratingFunctionTest, sum) {
+  GeneratingFunction gf1(2, {1, -1},
+                         {makeFracMatrix(2, 3, {{1, 2, 5}, {7, 2, 6}}),
+                          makeFracMatrix(2, 3, {{5, 2, 5}, {3, 7, 2}})},
+                         {{{3, 6}, {7, 2}}, {{2, 8}, {6, 3}}});
+  GeneratingFunction gf2(2, {1, 1},
+                         {makeFracMatrix(2, 3, {{6, 2, 1}, {4, 2, 6}}),
+                          makeFracMatrix(2, 3, {{3, 2, 6}, {9, 2, 5}})},
+                         {{{3, 7}, {5, 1}}, {{5, 2}, {6, 2}}});
+
+  GeneratingFunction sum = gf1 + gf2;
+  EXPECT_EQ_REPR_GENERATINGFUNCTION(
+      sum, GeneratingFunction(2, {1, -1, 1, 1},
+                              {makeFracMatrix(2, 3, {{1, 2, 5}, {7, 2, 6}}),
+                               makeFracMatrix(2, 3, {{5, 2, 5}, {3, 7, 2}}),
+                               makeFracMatrix(2, 3, {{6, 2, 1}, {4, 2, 6}}),
+                               makeFracMatrix(2, 3, {{3, 2, 6}, {9, 2, 5}})},
+                              {{{3, 6}, {7, 2}},
+                               {{2, 8}, {6, 3}},
+                               {{3, 7}, {5, 1}},
+                               {{5, 2}, {6, 2}}}));
+}
diff --git a/mlir/unittests/Analysis/Presburger/Utils.h b/mlir/unittests/Analysis/Presburger/Utils.h
index 2a9966c7ce2ea5..6b00898a7e2749 100644
--- a/mlir/unittests/Analysis/Presburger/Utils.h
+++ b/mlir/unittests/Analysis/Presburger/Utils.h
@@ -13,6 +13,7 @@
 #ifndef MLIR_UNITTESTS_ANALYSIS_PRESBURGER_UTILS_H
 #define MLIR_UNITTESTS_ANALYSIS_PRESBURGER_UTILS_H
 
+#include "mlir/Analysis/Presburger/GeneratingFunction.h"
 #include "mlir/Analysis/Presburger/IntegerRelation.h"
 #include "mlir/Analysis/Presburger/Matrix.h"
 #include "mlir/Analysis/Presburger/PWMAFunction.h"
@@ -72,9 +73,42 @@ inline void EXPECT_EQ_FRAC_MATRIX(FracMatrix a, FracMatrix b) {
       EXPECT_EQ(a(row, col), b(row, col));
 }
 
+// Check the coefficients (in order) of two generating functions.
+// Note that this is not a true equality check.
+inline void EXPECT_EQ_REPR_GENERATINGFUNCTION(detail::GeneratingFunction a,
+                                              detail::GeneratingFunction b) {
+  EXPECT_EQ(a.getNumParams(), b.getNumParams());
+
+  SmallVector<int> aSigns = a.getSigns();
+  SmallVector<int> bSigns = b.getSigns();
+  EXPECT_EQ(aSigns.size(), bSigns.size());
+  for (unsigned i = 0, e = aSigns.size(); i < e; i++)
+    EXPECT_EQ(aSigns[i], bSigns[i]);
+
+  std::vector<detail::ParamPoint> aNums = a.getNumerators();
+  std::vector<detail::ParamPoint> bNums = b.getNumerators();
+  EXPECT_EQ(aNums.size(), bNums.size());
+  for (unsigned i = 0, e = aNums.size(); i < e; i++)
+    EXPECT_EQ_FRAC_MATRIX(aNums[i], bNums[i]);
+
+  std::vector<std::vector<detail::Point>> aDens = a.getDenominators();
+  std::vector<std::vector<detail::Point>> bDens = b.getDenominators();
+  EXPECT_EQ(aDens.size(), bDens.size());
+  for (unsigned i = 0, e = aDens.size(); i < e; i++) {
+    EXPECT_EQ(aDens[i].size(), bDens[i].size());
+    for (unsigned j = 0, f = aDens[i].size(); j < f; j++) {
+      EXPECT_EQ(aDens[i][j].size(), bDens[i][j].size());
+      for (unsigned k = 0, g = aDens[i][j].size(); k < g; k++) {
+        EXPECT_EQ(aDens[i][j][k], bDens[i][j][k]);
+      }
+    }
+  }
+}
+
 // Check the coefficients (in order) of two quasipolynomials.
 // Note that this is not a true equality check.
-inline void EXPECT_EQ_REPR_QUASIPOLYNOMIAL(QuasiPolynomial a, QuasiPolynomial b) {
+inline void EXPECT_EQ_REPR_QUASIPOLYNOMIAL(QuasiPolynomial a,
+                                           QuasiPolynomial b) {
   EXPECT_EQ(a.getNumInputs(), b.getNumInputs());
 
   SmallVector<Fraction> aCoeffs = a.getCoefficients(),

``````````

</details>


https://github.com/llvm/llvm-project/pull/77199


More information about the cfe-commits mailing list