diff --git a/mlir/include/mlir/Analysis/Presburger/LinearTransform.h b/mlir/include/mlir/Analysis/Presburger/LinearTransform.h --- a/mlir/include/mlir/Analysis/Presburger/LinearTransform.h +++ b/mlir/include/mlir/Analysis/Presburger/LinearTransform.h @@ -32,7 +32,7 @@ // strictly below that of the previous column, and all columns which have only // zeros are at the end. static std::pair<unsigned, LinearTransform> - makeTransformToColumnEchelon(Matrix m); + makeTransformToColumnEchelon(const Matrix &m); // Returns an IntegerRelation having a constraint vector vT for every // constraint vector v in rel, where T is this transform. diff --git a/mlir/include/mlir/Analysis/Presburger/Matrix.h b/mlir/include/mlir/Analysis/Presburger/Matrix.h --- a/mlir/include/mlir/Analysis/Presburger/Matrix.h +++ b/mlir/include/mlir/Analysis/Presburger/Matrix.h @@ -154,6 +154,17 @@ SmallVector<int64_t, 8> postMultiplyWithColumn(ArrayRef<int64_t> colVec) const; + /// Given the current matrix M, returns the matrices H, U such that H is the + /// column hermite normal form of M, i.e. H = M * U, where U is unimodular and + /// the matrix H has the following restrictions: + /// - H is lower triangular. + /// - The leading coefficient (the first non-zero entry from the top, called + /// the pivot) of a non-zero column is always strictly below of the leading + /// coefficient of the column before it; moreover, it is positive. + /// - The elements to the right of pivots are zero and the elements to the + /// left of the pivots are nonnegative and strictly smaller than the pivot. + std::tuple<Matrix, Matrix> hermiteNormalForm() const; + /// Resize the matrix to the specified dimensions. If a dimension is smaller, /// the values are truncated; if it is bigger, the new values are initialized /// to zero. diff --git a/mlir/lib/Analysis/Presburger/LinearTransform.cpp b/mlir/lib/Analysis/Presburger/LinearTransform.cpp --- a/mlir/lib/Analysis/Presburger/LinearTransform.cpp +++ b/mlir/lib/Analysis/Presburger/LinearTransform.cpp @@ -15,101 +15,30 @@ LinearTransform::LinearTransform(Matrix &&oMatrix) : matrix(oMatrix) {} LinearTransform::LinearTransform(const Matrix &oMatrix) : matrix(oMatrix) {} -// Set M(row, targetCol) to its remainder on division by M(row, sourceCol) -// by subtracting from column targetCol an appropriate integer multiple of -// sourceCol. This brings M(row, targetCol) to the range [0, M(row, sourceCol)). -// Apply the same column operation to otherMatrix, with the same integer -// multiple. -static void modEntryColumnOperation(Matrix &m, unsigned row, unsigned sourceCol, - unsigned targetCol, Matrix &otherMatrix) { - assert(m(row, sourceCol) != 0 && "Cannot divide by zero!"); - assert((m(row, sourceCol) > 0 && m(row, targetCol) > 0) && - "Operands must be positive!"); - int64_t ratio = m(row, targetCol) / m(row, sourceCol); - m.addToColumn(sourceCol, targetCol, -ratio); - otherMatrix.addToColumn(sourceCol, targetCol, -ratio); -} - std::pair<unsigned, LinearTransform> -LinearTransform::makeTransformToColumnEchelon(Matrix m) { - // We start with an identity result matrix and perform operations on m - // until m is in column echelon form. We apply the same sequence of operations - // on resultMatrix to obtain a transform that takes m to column echelon - // form. - Matrix resultMatrix = Matrix::identity(m.getNumColumns()); - - unsigned echelonCol = 0; - // Invariant: in all rows above row, all columns from echelonCol onwards - // are all zero elements. In an iteration, if the curent row has any non-zero - // elements echelonCol onwards, we bring one to echelonCol and use it to - // make all elements echelonCol + 1 onwards zero. - for (unsigned row = 0; row < m.getNumRows(); ++row) { - // Search row for a non-empty entry, starting at echelonCol. - unsigned nonZeroCol = echelonCol; - for (unsigned e = m.getNumColumns(); nonZeroCol < e; ++nonZeroCol) { - if (m(row, nonZeroCol) == 0) - continue; - break; - } - - // Continue to the next row with the same echelonCol if this row is all - // zeros from echelonCol onwards. - if (nonZeroCol == m.getNumColumns()) - continue; - - // Bring the non-zero column to echelonCol. This doesn't affect rows - // above since they are all zero at these columns. - if (nonZeroCol != echelonCol) { - m.swapColumns(nonZeroCol, echelonCol); - resultMatrix.swapColumns(nonZeroCol, echelonCol); - } - - // Make m(row, echelonCol) non-negative. - if (m(row, echelonCol) < 0) { - m.negateColumn(echelonCol); - resultMatrix.negateColumn(echelonCol); - } - - // Make all the entries in row after echelonCol zero. - for (unsigned i = echelonCol + 1, e = m.getNumColumns(); i < e; ++i) { - // We make m(row, i) non-negative, and then apply the Euclidean GCD - // algorithm to (row, i) and (row, echelonCol). At the end, one of them - // has value equal to the gcd of the two entries, and the other is zero. - - if (m(row, i) < 0) { - m.negateColumn(i); - resultMatrix.negateColumn(i); - } - - unsigned targetCol = i, sourceCol = echelonCol; - // At every step, we set m(row, targetCol) %= m(row, sourceCol), and - // swap the indices sourceCol and targetCol. (not the columns themselves) - // This modulo is implemented as a subtraction - // m(row, targetCol) -= quotient * m(row, sourceCol), - // where quotient = floor(m(row, targetCol) / m(row, sourceCol)), - // which brings m(row, targetCol) to the range [0, m(row, sourceCol)). - // - // We are only allowed column operations; we perform the above - // for every row, i.e., the above subtraction is done as a column - // operation. This does not affect any rows above us since they are - // guaranteed to be zero at these columns. - while (m(row, targetCol) != 0 && m(row, sourceCol) != 0) { - modEntryColumnOperation(m, row, sourceCol, targetCol, resultMatrix); - std::swap(targetCol, sourceCol); - } - - // One of (row, echelonCol) and (row, i) is zero and the other is the gcd. - // Make it so that (row, echelonCol) holds the non-zero value. - if (m(row, echelonCol) == 0) { - m.swapColumns(i, echelonCol); - resultMatrix.swapColumns(i, echelonCol); +LinearTransform::makeTransformToColumnEchelon(const Matrix &m) { + // Compute the hermite normal form of m. The hermite form, is by defination + // column echelon. + auto [h, u] = m.hermiteNormalForm(); + + // Since the hermite normal form is lower triangular, a zero column means + // the rest of the columns are zero. Thus, once we find a zero column, we + // can stop. + unsigned col, e; + for (col = 0, e = m.getNumColumns(); col < e; ++col) { + bool zeroCol = true; + for (unsigned row = 0, f = m.getNumRows(); row < f; ++row) { + if (h(row, col) != 0) { + zeroCol = false; + break; } } - ++echelonCol; + if (zeroCol) + break; } - return {echelonCol, LinearTransform(std::move(resultMatrix))}; + return {col, LinearTransform(std::move(u))}; } IntegerRelation LinearTransform::applyTo(const IntegerRelation &rel) const { diff --git a/mlir/lib/Analysis/Presburger/Matrix.cpp b/mlir/lib/Analysis/Presburger/Matrix.cpp --- a/mlir/lib/Analysis/Presburger/Matrix.cpp +++ b/mlir/lib/Analysis/Presburger/Matrix.cpp @@ -247,6 +247,125 @@ return result; } +/// Set M(row, targetCol) to its remainder on division by M(row, sourceCol) +/// by subtracting from column targetCol an appropriate integer multiple of +/// sourceCol. This brings M(row, targetCol) to the range [0, M(row, +/// sourceCol)). Apply the same column operation to otherMatrix, with the same +/// integer multiple. +static void modEntryColumnOperation(Matrix &m, unsigned row, unsigned sourceCol, + unsigned targetCol, Matrix &otherMatrix) { + assert(m(row, sourceCol) != 0 && "Cannot divide by zero!"); + assert((m(row, sourceCol) > 0 && m(row, targetCol) > 0) && + "Operands must be positive!"); + int64_t ratio = m(row, targetCol) / m(row, sourceCol); + m.addToColumn(sourceCol, targetCol, -ratio); + otherMatrix.addToColumn(sourceCol, targetCol, -ratio); +} + +std::tuple<Matrix, Matrix> Matrix::hermiteNormalForm() const { + // We start with u as an identity matrix and perform operations on h until h + // is in hermite normal form. We apply the same sequence of operations on u to + // obtain a transform that takes h to hermite normal form. + Matrix h = *this; + Matrix u = Matrix::identity(h.getNumColumns()); + + unsigned echelonCol = 0; + // Invariant: in all rows above row, all columns from echelonCol onwards + // are all zero elements. In an iteration, if the curent row has any non-zero + // elements echelonCol onwards, we bring one to echelonCol and use it to + // make all elements echelonCol + 1 onwards zero. + for (unsigned row = 0; row < h.getNumRows(); ++row) { + // Search row for a non-empty entry, starting at echelonCol. + unsigned nonZeroCol = echelonCol; + for (unsigned e = h.getNumColumns(); nonZeroCol < e; ++nonZeroCol) { + if (h(row, nonZeroCol) == 0) + continue; + break; + } + + // Continue to the next row with the same echelonCol if this row is all + // zeros from echelonCol onwards. + if (nonZeroCol == h.getNumColumns()) + continue; + + // Bring the non-zero column to echelonCol. This doesn't affect rows + // above since they are all zero at these columns. + if (nonZeroCol != echelonCol) { + h.swapColumns(nonZeroCol, echelonCol); + u.swapColumns(nonZeroCol, echelonCol); + } + + // Make h(row, echelonCol) non-negative. + if (h(row, echelonCol) < 0) { + h.negateColumn(echelonCol); + u.negateColumn(echelonCol); + } + + // Make all the entries in row after echelonCol zero. + for (unsigned i = echelonCol + 1, e = h.getNumColumns(); i < e; ++i) { + // We make h(row, i) non-negative, and then apply the Euclidean GCD + // algorithm to (row, i) and (row, echelonCol). At the end, one of them + // has value equal to the gcd of the two entries, and the other is zero. + + if (h(row, i) < 0) { + h.negateColumn(i); + u.negateColumn(i); + } + + unsigned targetCol = i, sourceCol = echelonCol; + // At every step, we set h(row, targetCol) %= h(row, sourceCol), and + // swap the indices sourceCol and targetCol. (not the columns themselves) + // This modulo is implemented as a subtraction + // h(row, targetCol) -= quotient * h(row, sourceCol), + // where quotient = floor(h(row, targetCol) / h(row, sourceCol)), + // which brings h(row, targetCol) to the range [0, h(row, sourceCol)). + // + // We are only allowed column operations; we perform the above + // for every row, i.e., the above subtraction is done as a column + // operation. This does not affect any rows above us since they are + // guaranteed to be zero at these columns. + while (h(row, targetCol) != 0 && h(row, sourceCol) != 0) { + modEntryColumnOperation(h, row, sourceCol, targetCol, u); + std::swap(targetCol, sourceCol); + } + + // One of (row, echelonCol) and (row, i) is zero and the other is the gcd. + // Make it so that (row, echelonCol) holds the non-zero value. + if (h(row, echelonCol) == 0) { + h.swapColumns(i, echelonCol); + u.swapColumns(i, echelonCol); + } + } + + // Make all entries before echelonCol non-negative and strictly smaller + // than the pivot entry. + for (unsigned i = 0; i < echelonCol; ++i) { + // The scale to bring h(row, i) to the range [0, h(row, echelonCol)). + int64_t scale; + + if (h(row, i) < 0) { + // For negative entries, we add the smallest multiple of pivot that + // makes it non-negative. + int64_t absEntry = -h(row, i); + scale = absEntry / h(row, echelonCol); + if (absEntry % h(row, echelonCol) != 0) + scale++; + } else { + // For positive entries we subtract the largest multiple of pivot that + // makes it strictly smaller than pivot. + scale = -(h(row, i) / h(row, echelonCol)); + } + + h.addToColumn(echelonCol, i, scale); + u.addToColumn(echelonCol, i, scale); + } + + ++echelonCol; + } + + return {h, u}; +} + void Matrix::print(raw_ostream &os) const { for (unsigned row = 0; row < nRows; ++row) { for (unsigned column = 0; column < nColumns; ++column) diff --git a/mlir/unittests/Analysis/Presburger/MatrixTest.cpp b/mlir/unittests/Analysis/Presburger/MatrixTest.cpp --- a/mlir/unittests/Analysis/Presburger/MatrixTest.cpp +++ b/mlir/unittests/Analysis/Presburger/MatrixTest.cpp @@ -7,6 +7,7 @@ //===----------------------------------------------------------------------===// #include "mlir/Analysis/Presburger/Matrix.h" +#include "./Utils.h" #include <gmock/gmock.h> #include <gtest/gtest.h> @@ -191,3 +192,78 @@ for (unsigned col = 0; col < 7; ++col) EXPECT_EQ(mat(row, col), row >= 3 || col >= 3 ? 0 : int(10 * row + col)); } + +TEST(MatrixTest, hermiteNormalForm) { + // TODO: Add a check to test the original statement of hermite normal form + // instead of using a precomputed result. + + { + // Hermite form of a unimodular is the identity matrix. + Matrix mat = makeMatrix(3, 3, {{2, 3, 6}, {3, 2, 3}, {17, 11, 16}}); + Matrix hermiteForm = makeMatrix(3, 3, {{1, 0, 0}, {0, 1, 0}, {0, 0, 1}}); + + auto [h, u] = mat.hermiteNormalForm(); + + for (unsigned row = 0; row < 3; row++) + for (unsigned col = 0; col < 3; col++) + EXPECT_EQ(h(row, col), hermiteForm(row, col)); + } + + { + // Hermite form of a unimodular is the identity matrix. + Matrix mat = makeMatrix( + 4, 4, + {{-6, -1, -19, -20}, {0, 1, 0, 0}, {-5, 0, -15, -16}, {6, 0, 18, 19}}); + + Matrix hermiteForm = makeMatrix( + 4, 4, {{1, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, 1, 0}, {0, 0, 0, 1}}); + + auto [h, u] = mat.hermiteNormalForm(); + + for (unsigned row = 0; row < 4; row++) + for (unsigned col = 0; col < 4; col++) + EXPECT_EQ(h(row, col), hermiteForm(row, col)); + } + + { + Matrix mat = makeMatrix( + 4, 4, {{3, 3, 1, 4}, {0, 1, 0, 0}, {0, 0, 19, 16}, {0, 0, 0, 3}}); + + Matrix hermiteForm = makeMatrix( + 4, 4, {{1, 0, 0, 0}, {0, 1, 0, 0}, {1, 0, 3, 0}, {18, 0, 54, 57}}); + + auto [h, u] = mat.hermiteNormalForm(); + + for (unsigned row = 0; row < 4; row++) + for (unsigned col = 0; col < 4; col++) + EXPECT_EQ(h(row, col), hermiteForm(row, col)); + } + + { + Matrix mat = makeMatrix( + 4, 4, {{3, 3, 1, 4}, {0, 1, 0, 0}, {0, 0, 19, 16}, {0, 0, 0, 3}}); + + Matrix hermiteForm = makeMatrix( + 4, 4, {{1, 0, 0, 0}, {0, 1, 0, 0}, {1, 0, 3, 0}, {18, 0, 54, 57}}); + + auto [h, u] = mat.hermiteNormalForm(); + + for (unsigned row = 0; row < 4; row++) + for (unsigned col = 0; col < 4; col++) + EXPECT_EQ(h(row, col), hermiteForm(row, col)); + } + + { + Matrix mat = + makeMatrix(3, 5, {{0, 2, 0, 7, 1}, {-1, 0, 0, -3, 0}, {0, 4, 1, 0, 8}}); + + Matrix hermiteForm = + makeMatrix(3, 5, {{1, 0, 0, 0, 0}, {0, 1, 0, 0, 0}, {0, 0, 1, 0, 0}}); + + auto [h, u] = mat.hermiteNormalForm(); + + for (unsigned row = 0; row < 3; row++) + for (unsigned col = 0; col < 5; col++) + EXPECT_EQ(h(row, col), hermiteForm(row, col)); + } +}