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));
+  }
+}