[Mlir-commits] [mlir] [MLIR][Presburger] Add support for Smith normal form (PR #185328)
llvmlistbot at llvm.org
llvmlistbot at llvm.org
Sun Mar 8 15:38:44 PDT 2026
llvmbot wrote:
<!--LLVM PR SUMMARY COMMENT-->
@llvm/pr-subscribers-mlir-presburger
@llvm/pr-subscribers-mlir
Author: Yue Huang (AdUhTkJm)
<details>
<summary>Changes</summary>
FPL already has support for computing Hermite normal form for integer matrices. Here we add support to computing Smith normal form.
This is a preparation for Barvinok's algorithm, where we must find particular solution in order to project lower-dimensional polyhedra into full-dimensional ones.
The implementation here follows the algorithm in [wikipedia](https://en.wikipedia.org/wiki/Smith_normal_form#Algorithm).
---
Full diff: https://github.com/llvm/llvm-project/pull/185328.diff
3 Files Affected:
- (modified) mlir/include/mlir/Analysis/Presburger/Matrix.h (+9)
- (modified) mlir/lib/Analysis/Presburger/Matrix.cpp (+114)
- (modified) mlir/unittests/Analysis/Presburger/MatrixTest.cpp (+50)
``````````diff
diff --git a/mlir/include/mlir/Analysis/Presburger/Matrix.h b/mlir/include/mlir/Analysis/Presburger/Matrix.h
index 15069fa2a55f1..7c24b9c8d4c40 100644
--- a/mlir/include/mlir/Analysis/Presburger/Matrix.h
+++ b/mlir/include/mlir/Analysis/Presburger/Matrix.h
@@ -274,6 +274,15 @@ class IntMatrix : public Matrix<DynamicAPInt> {
/// pivot.
std::pair<IntMatrix, IntMatrix> computeHermiteNormalForm() const;
+ /// Given the current matrix M, returns the matrices U, D, V such that
+ /// UMV = D, where D is called the Smith Normal Form (SNF).
+ /// The matrices have the following properties:
+ /// - U, V are unimodular. In other words, det(U), det(V) are 1 or -1;
+ /// their inverse also contain integer entries.
+ /// - D is diagonal.
+ /// - For all i, the diagonal element D_{i, i} divides D_{i + 1, i + 1}.
+ std::tuple<IntMatrix, IntMatrix, IntMatrix> computeSmithNormalForm() const;
+
/// Divide the first `nCols` of the specified row by their GCD.
/// Returns the GCD of the first `nCols` of the specified row.
DynamicAPInt normalizeRow(unsigned row, unsigned nCols);
diff --git a/mlir/lib/Analysis/Presburger/Matrix.cpp b/mlir/lib/Analysis/Presburger/Matrix.cpp
index 83a2c280c3d4e..4b178d0ccdacf 100644
--- a/mlir/lib/Analysis/Presburger/Matrix.cpp
+++ b/mlir/lib/Analysis/Presburger/Matrix.cpp
@@ -535,6 +535,120 @@ std::pair<IntMatrix, IntMatrix> IntMatrix::computeHermiteNormalForm() const {
return {h, u};
}
+std::tuple<IntMatrix, IntMatrix, IntMatrix>
+IntMatrix::computeSmithNormalForm() const {
+ IntMatrix d = *this;
+ // The matrix U records row operations, while V records column operations.
+ // When we operate on the matrix D in the following code, we record the same
+ // operations on U and V respectively.
+ IntMatrix u = IntMatrix::identity(d.getNumRows());
+ IntMatrix v = IntMatrix::identity(d.getNumColumns());
+
+ unsigned numRows = d.getNumRows();
+ unsigned numCols = d.getNumColumns();
+ for (unsigned i = 0, e = std::min(numRows, numCols); i < e; i++) {
+ bool changed;
+ do {
+ changed = false;
+
+ // Find the entry in the submatrix d(i:, i:) with the smallest non-zero
+ // absolute value.
+ // The element is denoted d(p, q), and is the pivot.
+ unsigned p = i, q = i;
+ DynamicAPInt minVal(0);
+ for (unsigned r = i; r < numRows; r++)
+ for (unsigned c = i; c < numCols; c++) {
+ DynamicAPInt val = llvm::abs(d(r, c));
+ if (val != 0 && (minVal == 0 || val < minVal)) {
+ minVal = val;
+ p = r;
+ q = c;
+ }
+ }
+
+ // The remaining submatrix is zero.
+ if (minVal == 0)
+ break;
+
+ // Bring pivot to d(i, i). Record the operation in u, v respectively.
+ if (p != i) {
+ d.swapRows(p, i);
+ u.swapRows(p, i);
+ }
+ if (q != i) {
+ d.swapColumns(q, i);
+ v.swapColumns(q, i);
+ }
+
+ // The pivot should be positive.
+ if (d(i, i) < 0) {
+ d.negateRow(i);
+ u.negateRow(i);
+ }
+
+ DynamicAPInt pivot = d(i, i);
+
+ // Clear other entries in row i and column i with Euclid's algorithm.
+ for (unsigned r = 0; r < numRows; ++r) {
+ if (r == i)
+ continue;
+ while (d(r, i) != 0) {
+ auto quotient = d(r, i) / d(i, i);
+ d.addToRow(i, r, -quotient);
+ u.addToRow(i, r, -quotient);
+
+ if (llvm::abs(d(r, i)) < llvm::abs(d(i, i)) && d(r, i) != 0) {
+ d.swapRows(r, i);
+ u.swapRows(r, i);
+ changed = true;
+ }
+ }
+ }
+ // Similar to the rows operations, this time it works on columns.
+ for (unsigned c = 0; c < numCols; ++c) {
+ if (c == i)
+ continue;
+ while (d(i, c) != 0) {
+ DynamicAPInt quotient = d(i, c) / d(i, i);
+
+ d.addToColumn(i, c, -quotient);
+ v.addToColumn(i, c, -quotient);
+
+ if (llvm::abs(d(i, c)) < llvm::abs(d(i, i)) && d(i, c) != 0) {
+ d.swapColumns(c, i);
+ v.swapColumns(c, i);
+ changed = true;
+ }
+ assert(d(i, i) != 0);
+ }
+ }
+
+ for (unsigned r = i + 1; r < numRows; ++r) {
+ bool breakFlag = false;
+ for (unsigned c = i + 1; c < numCols; ++c) {
+ if (d(r, c) % pivot != 0) {
+ // Add row r to row i. This brings d(r, c) into the i-th row,
+ // creating a new value at d(i, c) that will be used to reduce the
+ // pivot size.
+ d.addToRow(r, i, 1);
+ u.addToRow(r, i, 1);
+ changed = true;
+
+ // We must break and restart the inner process to find the new,
+ // smaller pivot in the submatrix d(i:, i:).
+ breakFlag = true;
+ break;
+ }
+ }
+ if (breakFlag)
+ break;
+ }
+ } while (changed);
+ }
+
+ return {u, d, v};
+}
+
DynamicAPInt IntMatrix::normalizeRow(unsigned row, unsigned cols) {
return normalizeRange(getRow(row).slice(0, cols));
}
diff --git a/mlir/unittests/Analysis/Presburger/MatrixTest.cpp b/mlir/unittests/Analysis/Presburger/MatrixTest.cpp
index e2c2a9bcb7d26..db4b502affcaf 100644
--- a/mlir/unittests/Analysis/Presburger/MatrixTest.cpp
+++ b/mlir/unittests/Analysis/Presburger/MatrixTest.cpp
@@ -258,6 +258,56 @@ TEST(MatrixTest, computeHermiteNormalForm) {
}
}
+static IntMatrix matmul(const IntMatrix &a, const IntMatrix &b) {
+ assert(a.getNumColumns() == b.getNumRows());
+ unsigned n = a.getNumRows();
+ unsigned m = b.getNumRows();
+ unsigned p = b.getNumColumns();
+ IntMatrix result(n, p);
+
+ for (unsigned i = 0; i < n; i++) {
+ for (unsigned j = 0; j < m; j++) {
+ for (unsigned k = 0; k < p; k++) {
+ result(i, k) += a(i, j) * b(j, k);
+ }
+ }
+ }
+ return result;
+}
+
+TEST(MatrixTest, computeSmithNormalForm) {
+ IntMatrix mat =
+ makeIntMatrix(4, 3, {{2, 5, 8}, {3, 6, 9}, {4, 7, 1}, {5, 8, 2}});
+
+ auto [u, d, v] = mat.computeSmithNormalForm();
+
+ // Check u and v are unimodular.
+ EXPECT_TRUE(llvm::abs(u.determinant()) == 1);
+ EXPECT_TRUE(llvm::abs(v.determinant()) == 1);
+
+ // Check u @ mat @ v = d (@ is matrix multiplication).
+ EXPECT_TRUE(matmul(u, matmul(mat, v)) == d);
+
+ // Check d is diagonal.
+ for (unsigned i = 0, e = d.getNumRows(); i < e; i++) {
+ for (unsigned j = 0, f = d.getNumColumns(); j < f; j++) {
+ if (i == j)
+ EXPECT_TRUE(d(i, j) != 0);
+ else
+ EXPECT_TRUE(d(i, j) == 0);
+ }
+ }
+
+ // Check d(i, i) divides d(i + 1, i + 1).
+ unsigned end = std::min(d.getNumRows(), d.getNumColumns()) - 1;
+ for (unsigned i = 0; i < end; i++) {
+ if (d(i, i) == 0)
+ break;
+
+ EXPECT_TRUE(d(i + 1, i + 1) % d(i, i) == 0);
+ }
+}
+
TEST(MatrixTest, inverse) {
IntMatrix mat1 = makeIntMatrix(2, 2, {{2, 1}, {7, 0}});
EXPECT_EQ(mat1.determinant(), -7);
``````````
</details>
https://github.com/llvm/llvm-project/pull/185328
More information about the Mlir-commits
mailing list