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 - 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 @@ -156,6 +156,18 @@ SmallVector postMultiplyWithColumn(ArrayRef 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 the pivots are zero and the elements to + /// the left of the pivots are nonnegative and strictly smaller than the + /// pivot. + std::pair computeHermiteNormalForm() 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,29 @@ 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 -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. This, is by definition, is in column + // echelon form. + auto [h, u] = m.computeHermiteNormalForm(); + + // Since the matrix is in column ecehlon form, 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 @@ -251,6 +251,106 @@ 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 && "Source must be positive!"); + int64_t ratio = -floorDiv(m(row, targetCol), m(row, sourceCol)); + m.addToColumn(sourceCol, targetCol, ratio); + otherMatrix.addToColumn(sourceCol, targetCol, ratio); +} + +std::pair Matrix::computeHermiteNormalForm() 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) + modEntryColumnOperation(h, row, echelonCol, i, u); + + ++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 #include @@ -191,3 +192,58 @@ for (unsigned col = 0; col < 7; ++col) EXPECT_EQ(mat(row, col), row >= 3 || col >= 3 ? 0 : int(10 * row + col)); } + +static void checkHermiteNormalForm(const Matrix &mat, + const Matrix &hermiteForm) { + auto [h, u] = mat.computeHermiteNormalForm(); + + for (unsigned row = 0; row < mat.getNumRows(); row++) + for (unsigned col = 0; col < mat.getNumColumns(); col++) + EXPECT_EQ(h(row, col), hermiteForm(row, col)); +} + +TEST(MatrixTest, computeHermiteNormalForm) { + // TODO: Add a check to test the original statement of hermite normal form + // instead of using a precomputed result. + + { + // Hermite form of a unimodular matrix 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}}); + checkHermiteNormalForm(mat, hermiteForm); + } + + { + // 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}}); + checkHermiteNormalForm(mat, hermiteForm); + } + + { + 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}}); + checkHermiteNormalForm(mat, hermiteForm); + } + + { + 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}}); + checkHermiteNormalForm(mat, hermiteForm); + } + + { + 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}}); + checkHermiteNormalForm(mat, hermiteForm); + } +}