Template Matrix to allow both MPInt and Fraction

[Abhinav271828]

Summary

Template Matrix to Matrix<T> (for MPInt and Fraction) with explicit instantiation
Duplicate makeMatrix to makeIntMatrix and makeFracMatrix
Implement arithmetic operations for Fraction for compatibility

[Abhinav271828]

Template Matrix to work with MPInt and Fraction
Add printing capability to Fraction
All previous uses of Matrix are replaced with Matrix<MPInt>

[Groverkss]

Looks like you forgot some commits. Could you check?

[Abhinav271828]

In D156906#4600783, @Groverkss wrote:

Looks like you forgot some commits. Could you check?

Changed the diff, does it look all right now?

[arjunp]

I suggest breaking up this patch further. Let's keep this patch to only templating the Matrix class. Adding new functionality can be done in a separate patch. How does inverse work when it's a matrix of integers?

[Abhinav271828]

In D156906#4624362, @arjunp wrote:

I suggest breaking up this patch further. Let's keep this patch to only templating the Matrix class. Adding new functionality can be done in a separate patch. How does inverse work when it's a matrix of integers?

Done! Only the templating is part of the patch now.
Inverse is only defined for rational matrices as of now.

[arjunp]

Templated C++ classes usually need to have all their implementations in the header but this increases compile time. However we only need Matrix<MPInt> and Matrix<Fraction>.

I think the best option here is to just rename Matrix<MPInt> to IntMatrix (using IntMatrix = detail::Matrix<MPInt>), hide the templated matrix class in a namespace detail, and document that it's not intended for general use. You then have to explicitly instantiate Matrix<MPInt> in the .cpp. (forget about Matrix<Fraction> for now.) Add this line to the end of the .cpp

template class Matix<MPInt>;

(https://stackoverflow.com/questions/2351148/explicit-template-instantiation-when-is-it-used)

(I also considered putting static_assert that the type is a supported one, but this seems cleaner to me.)

In D156906#4624875, @Abhinav271828 wrote:

In D156906#4624362, @arjunp wrote:

I suggest breaking up this patch further. Let's keep this patch to only templating the Matrix class. Adding new functionality can be done in a separate patch. How does inverse work when it's a matrix of integers?

Done! Only the templating is part of the patch now.
Inverse is only defined for rational matrices as of now.

It's not a good idea to list functions in the header that then surprisingly turn out to not have been implemented. I would probably just move it out into a free function for that patch. The Matrix class atm is supposed to be a 2D array anyway, not a linear transform class. We can discuss whether that's a good idea but it's beyond the scope of these patches.

mlir/lib/Analysis/Presburger/Matrix.cpp

240

why was this changed?

252

same here

mlir/unittests/Analysis/Presburger/MatrixTest.cpp

212–213

Please rename makeMatrix to makeIntMatrix and use that instead of templating makeMatrix and adding constructors hre. The current solution is difficult to read

[Abhinav271828]

mlir/lib/Analysis/Presburger/Matrix.cpp

240

Fraction doesn't have += implemented, so the compilation failed for that instantiation.

mlir/unittests/Analysis/Presburger/MatrixTest.cpp

212–213

What about when we need to use makeMatrix for fractions? Should I duplicated makeMatrix to makeIntMatrix and makeFracMatrix?

[arjunp]

mlir/lib/Analysis/Presburger/Matrix.cpp

380–387

The braces should be in the same line and namespace contents shouldn't be indented

mlir/unittests/Analysis/Presburger/Utils.h

43–45

Not sure if this will work, you might need an std::pair to make the brace initialization actually work

[arjunp]

mlir/include/mlir/Analysis/Presburger/Matrix.h

33

Please mention here that this is only intended for use with Fraction and MPInt. And add a static assert to the top of the class checking that T is one of these types.

static_assert(std::is_same_v<T,MPInt> || std::is_same_v<T,Fraction>, "T must be MPInt or Fraction.")

[arjunp]

Please write full sentences in the patch summary.

[Abhinav271828]

mlir/unittests/Analysis/Presburger/Utils.h

43–45

The tests use the constructor itself, e.g. makeFracMatrix(2, 2, {{Fraction(1, 1), Fraction(0, 1)}, ...})

[Groverkss]

This patch was sent as a pull request instead: https://github.com/llvm/llvm-project/pull/65272. Please abandon this revision to reduce the number of open revisions we have.

mlir/include/mlir/Analysis/Presburger/Fraction.h

 #define MLIR_ANALYSIS_PRESBURGER_FRACTION_H
 
 #include "mlir/Analysis/Presburger/MPInt.h"
+#include "mlir/Analysis/Presburger/Utils.h"
 #include "mlir/Support/MathExtras.h"
 
 namespace mlir {

   Fraction() = default;
 
   /// Construct a Fraction from a numerator and denominator.
-  Fraction(const MPInt &oNum, const MPInt &oDen) : num(oNum), den(oDen) {
+  Fraction(const MPInt &oNum, const MPInt &oDen = MPInt(1)) : num(oNum), den(oDen) {
     if (den < 0) {
       num = -num;
       den = -den;
     }
   }
   /// Overloads for passing literals.
-  Fraction(const MPInt &num, int64_t den) : Fraction(num, MPInt(den)) {}
-  Fraction(int64_t num, const MPInt &den) : Fraction(MPInt(num), den) {}
+  Fraction(const MPInt &num, int64_t den = 1) : Fraction(num, MPInt(den)) {}
+  Fraction(int64_t num, const MPInt &den = MPInt(1)) : Fraction(MPInt(num), den) {}
   Fraction(int64_t num, int64_t den) : Fraction(MPInt(num), MPInt(den)) {}
 
   // Return the value of the fraction as an integer. This should only be called

     return num / den;
   }
 
+  llvm::raw_ostream &print(llvm::raw_ostream &os) const {
+    return os << "(" << num << "/" << den << ")";
+  }
+
   /// The numerator and denominator, respectively. The denominator is always
   /// positive.
   MPInt num{0}, den{1};

   return Fraction(x.num * y.num, x.den * y.den);
 }
 
+inline Fraction operator/(const Fraction &x, const Fraction &y) {
+  return Fraction(x.num * y.den, x.den * y.num);
+}
+
+inline Fraction operator+(const Fraction &x, const Fraction &y) {
+  return Fraction(x.num * y.den + x.den * y.num, x.den * y.den);
+}
+
+inline Fraction operator-(const Fraction &x, const Fraction &y) {
+  return Fraction(x.num * y.den - x.den * y.num, x.den * y.den);
+}
+
+inline llvm::raw_ostream &operator<<(llvm::raw_ostream &os, const Fraction &x) {
+  x.print(os);
+  return os;
+}
+
 } // namespace presburger
 } // namespace mlir
 

mlir/include/mlir/Analysis/Presburger/IntegerRelation.h

   /// bounded. The span of the returned vectors is guaranteed to contain all
   /// such vectors. The returned vectors are NOT guaranteed to be linearly
   /// independent. This function should not be called on empty sets.
-  Matrix getBoundedDirections() const;
+  Matrix<MPInt> getBoundedDirections() const;
 
   /// Find an integer sample point satisfying the constraints using a
   /// branch and bound algorithm with generalized basis reduction, with some

   PresburgerSpace space;
 
   /// Coefficients of affine equalities (in == 0 form).
-  Matrix equalities;
+  Matrix<MPInt> equalities;
 
   /// Coefficients of affine inequalities (in >= 0 form).
-  Matrix inequalities;
+  Matrix<MPInt> inequalities;
 };
 
 /// An IntegerPolyhedron represents the set of points from a PresburgerSpace

mlir/include/mlir/Analysis/Presburger/LinearTransform.h

 
 class LinearTransform {
 public:
-  explicit LinearTransform(Matrix &&oMatrix);
-  explicit LinearTransform(const Matrix &oMatrix);
+  explicit LinearTransform(Matrix<MPInt> &&oMatrix);
+  explicit LinearTransform(const Matrix<MPInt> &oMatrix);
 
   // Returns a linear transform T such that MT is M in column echelon form.
   // Also returns the number of non-zero columns in MT.

   // strictly below that of the previous column, and all columns which have only
   // zeros are at the end.
   static std::pair<unsigned, LinearTransform>
-  makeTransformToColumnEchelon(const Matrix &m);
+  makeTransformToColumnEchelon(const Matrix<MPInt> &m);
 
   // Returns an IntegerRelation having a constraint vector vT for every
   // constraint vector v in rel, where T is this transform.

     return matrix.postMultiplyWithColumn(colVec);
   }
 
+  // Compute the determinant of the transform by converting it to row echelon
+  // form and then taking the product of the diagonal.
+  MPInt determinant();
+
 private:
-  Matrix matrix;
+  Matrix<MPInt> matrix;
 };
 
 } // namespace presburger

mlir/include/mlir/Analysis/Presburger/Matrix.h

 #ifndef MLIR_ANALYSIS_PRESBURGER_MATRIX_H
 #define MLIR_ANALYSIS_PRESBURGER_MATRIX_H
 
-#include "mlir/Analysis/Presburger/MPInt.h"
 #include "mlir/Support/LLVM.h"
 #include "llvm/ADT/ArrayRef.h"
 #include "llvm/Support/raw_ostream.h"

 /// (i, j) is stored at data[i*nReservedColumns + j]. The reserved but unused
 /// columns always have all zero values. The reserved rows are just reserved
 /// space in the underlying SmallVector's capacity.
+template<typename T>
 class Matrix {
 public:
   Matrix() = delete;

   static Matrix identity(unsigned dimension);
 
   /// Access the element at the specified row and column.
-  MPInt &at(unsigned row, unsigned column) {
+  T &at(unsigned row, unsigned column) {
     assert(row < nRows && "Row outside of range");
     assert(column < nColumns && "Column outside of range");
     return data[row * nReservedColumns + column];
   }
 
-  MPInt at(unsigned row, unsigned column) const {
+  T at(unsigned row, unsigned column) const {
     assert(row < nRows && "Row outside of range");
     assert(column < nColumns && "Column outside of range");
     return data[row * nReservedColumns + column];
   }
 
-  MPInt &operator()(unsigned row, unsigned column) { return at(row, column); }
+  T &operator()(unsigned row, unsigned column) { return at(row, column); }
 
-  MPInt operator()(unsigned row, unsigned column) const {
+  T operator()(unsigned row, unsigned column) const {
     return at(row, column);
   }
 

   void reserveRows(unsigned rows);
 
   /// Get a [Mutable]ArrayRef corresponding to the specified row.
-  MutableArrayRef<MPInt> getRow(unsigned row);
-  ArrayRef<MPInt> getRow(unsigned row) const;
+  MutableArrayRef<T> getRow(unsigned row);
+  ArrayRef<T> getRow(unsigned row) const;
 
   /// Set the specified row to `elems`.
-  void setRow(unsigned row, ArrayRef<MPInt> elems);
+  void setRow(unsigned row, ArrayRef<T> elems);
 
   /// Insert columns having positions pos, pos + 1, ... pos + count - 1.
   /// Columns that were at positions 0 to pos - 1 will stay where they are;

 
   void copyRow(unsigned sourceRow, unsigned targetRow);
 
-  void fillRow(unsigned row, const MPInt &value);
-  void fillRow(unsigned row, int64_t value) { fillRow(row, MPInt(value)); }
+  void fillRow(unsigned row, const T &value);
+  void fillRow(unsigned row, int64_t value) { fillRow(row, T(value)); }
 
   /// Add `scale` multiples of the source row to the target row.
-  void addToRow(unsigned sourceRow, unsigned targetRow, const MPInt &scale);
+  void addToRow(unsigned sourceRow, unsigned targetRow, const T &scale);
   void addToRow(unsigned sourceRow, unsigned targetRow, int64_t scale) {
-    addToRow(sourceRow, targetRow, MPInt(scale));
+    addToRow(sourceRow, targetRow, T(scale));
   }
   /// Add `scale` multiples of the rowVec row to the specified row.
-  void addToRow(unsigned row, ArrayRef<MPInt> rowVec, const MPInt &scale);
+  void addToRow(unsigned row, ArrayRef<T> rowVec, const T &scale);
 
   /// Add `scale` multiples of the source column to the target column.
   void addToColumn(unsigned sourceColumn, unsigned targetColumn,
-                   const MPInt &scale);
+                   const T &scale);
   void addToColumn(unsigned sourceColumn, unsigned targetColumn,
                    int64_t scale) {
-    addToColumn(sourceColumn, targetColumn, MPInt(scale));
+    addToColumn(sourceColumn, targetColumn, T(scale));
   }
 
   /// Negate the specified column.

 
   /// Divide the first `nCols` of the specified row by their GCD.
   /// Returns the GCD of the first `nCols` of the specified row.
-  MPInt normalizeRow(unsigned row, unsigned nCols);
+  T normalizeRow(unsigned row, unsigned nCols);
   /// Divide the columns of the specified row by their GCD.
   /// Returns the GCD of the columns of the specified row.
-  MPInt normalizeRow(unsigned row);
+  T normalizeRow(unsigned row);
 
   /// The given vector is interpreted as a row vector v. Post-multiply v with
   /// this matrix, say M, and return vM.
-  SmallVector<MPInt, 8> preMultiplyWithRow(ArrayRef<MPInt> rowVec) const;
+  SmallVector<T, 8> preMultiplyWithRow(ArrayRef<T> rowVec) const;
 
   /// The given vector is interpreted as a column vector v. Pre-multiply v with
   /// this matrix, say M, and return Mv.
-  SmallVector<MPInt, 8> postMultiplyWithColumn(ArrayRef<MPInt> colVec) const;
+  SmallVector<T, 8> postMultiplyWithColumn(ArrayRef<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

   unsigned appendExtraRow();
   /// Same as above, but copy the given elements into the row. The length of
   /// `elems` must be equal to the number of columns.
-  unsigned appendExtraRow(ArrayRef<MPInt> elems);
+  unsigned appendExtraRow(ArrayRef<T> elems);
 
   /// Print the matrix.
   void print(raw_ostream &os) const;

 
   /// Stores the data. data.size() is equal to nRows * nReservedColumns.
   /// data.capacity() / nReservedColumns is the number of reserved rows.
-  SmallVector<MPInt, 16> data;
+  SmallVector<T, 16> data;
 };
 
 } // namespace presburger

mlir/include/mlir/Analysis/Presburger/PWMAFunction.h

 /// value of the function at a specified point.
 class MultiAffineFunction {
 public:
-  MultiAffineFunction(const PresburgerSpace &space, const Matrix &output)
+  MultiAffineFunction(const PresburgerSpace &space, const Matrix<MPInt> &output)
       : space(space), output(output),
         divs(space.getNumVars() - space.getNumRangeVars()) {
     assertIsConsistent();
   }
 
-  MultiAffineFunction(const PresburgerSpace &space, const Matrix &output,
+  MultiAffineFunction(const PresburgerSpace &space, const Matrix<MPInt> &output,
                       const DivisionRepr &divs)
       : space(space), output(output), divs(divs) {
     assertIsConsistent();

   PresburgerSpace getOutputSpace() const { return space.getRangeSpace(); }
 
   /// Get a matrix with each row representing row^th output expression.
-  const Matrix &getOutputMatrix() const { return output; }
+  const Matrix<MPInt> &getOutputMatrix() const { return output; }
   /// Get the `i^th` output expression.
   ArrayRef<MPInt> getOutputExpr(unsigned i) const { return output.getRow(i); }
 

   /// The function's output is a tuple of integers, with the ith element of the
   /// tuple defined by the affine expression given by the ith row of this output
   /// matrix.
-  Matrix output;
+  Matrix<MPInt> output;
 
   /// Storage for division representation for each local variable in space.
   DivisionRepr divs;

mlir/include/mlir/Analysis/Presburger/Simplex.h

   unsigned nSymbol;
 
   /// The matrix representing the tableau.
-  Matrix tableau;
+  Matrix<MPInt> tableau;
 
   /// This is true if the tableau has been detected to be empty, false
   /// otherwise.

 
   /// Reduce the given basis, starting at the specified level, using general
   /// basis reduction.
-  void reduceBasis(Matrix &basis, unsigned level);
+  void reduceBasis(Matrix<MPInt> &basis, unsigned level);
 };
 
 /// Takes a snapshot of the simplex state on construction and rolls back to the

mlir/include/mlir/Analysis/Presburger/Utils.h

   /// Each row of the Matrix represents a single division dividend. The
   /// `i^th` row represents the dividend of the variable at `divOffset + i`
   /// in the constraint system (and the `i^th` local variable).
-  Matrix dividends;
+  Matrix<MPInt> dividends;
 
   /// Denominators of each division. If a denominator of a division is `0`, the
   /// division variable is considered to not have a division representation.

mlir/lib/Analysis/FlatLinearValueConstraints.cpp

          "AffineMap cannot produce divs without local representation");
 
   // TODO: We shouldn't have to do this conversion.
-  Matrix mat(map.getNumResults(), map.getNumInputs() + divs.getNumDivs() + 1);
+  Matrix<MPInt> mat(map.getNumResults(), map.getNumInputs() + divs.getNumDivs() + 1);
   for (unsigned i = 0, e = flattenedExprs.size(); i < e; ++i)
     for (unsigned j = 0, f = flattenedExprs[i].size(); j < f; ++j)
       mat(i, j) = flattenedExprs[i][j];

mlir/lib/Analysis/Presburger/IntegerRelation.cpp

 //
 // It is sufficient to check the perpendiculars of the constraints, as the set
 // of perpendiculars which are bounded must span all bounded directions.
-Matrix IntegerRelation::getBoundedDirections() const {
+Matrix<MPInt> IntegerRelation::getBoundedDirections() const {
   // Note that it is necessary to add the equalities too (which the constructor
   // does) even though we don't need to check if they are bounded; whether an
   // inequality is bounded or not depends on what other constraints, including

   // The direction vector is given by the coefficients and does not include the
   // constant term, so the matrix has one fewer column.
   unsigned dirsNumCols = getNumCols() - 1;
-  Matrix dirs(boundedIneqs.size() + getNumEqualities(), dirsNumCols);
+  Matrix<MPInt> dirs(boundedIneqs.size() + getNumEqualities(), dirsNumCols);
 
   // Copy the bounded inequalities.
   unsigned row = 0;

   // m is a matrix containing, in each row, a vector in which S is
   // bounded, such that the linear span of all these dimensions contains all
   // bounded dimensions in S.
-  Matrix m = getBoundedDirections();
+  Matrix<MPInt> m = getBoundedDirections();
   // In column echelon form, each row of m occupies only the first rank(m)
   // columns and has zeros on the other columns. The transform T that brings S
   // to column echelon form is unimodular as well, so this is a suitable

mlir/lib/Analysis/Presburger/LinearTransform.cpp

 using namespace mlir;
 using namespace presburger;
 
-LinearTransform::LinearTransform(Matrix &&oMatrix) : matrix(oMatrix) {}
-LinearTransform::LinearTransform(const Matrix &oMatrix) : matrix(oMatrix) {}
+LinearTransform::LinearTransform(Matrix<MPInt> &&oMatrix) : matrix(oMatrix) {}
+LinearTransform::LinearTransform(const Matrix<MPInt> &oMatrix) : matrix(oMatrix) {}
 
 std::pair<unsigned, LinearTransform>
-LinearTransform::makeTransformToColumnEchelon(const Matrix &m) {
+LinearTransform::makeTransformToColumnEchelon(const Matrix<MPInt> &m) {
   // Compute the hermite normal form of m. This, is by definition, is in column
   // echelon form.
   auto [h, u] = m.computeHermiteNormalForm();

mlir/lib/Analysis/Presburger/Matrix.cpp

 //===----------------------------------------------------------------------===//
 
 #include "mlir/Analysis/Presburger/Matrix.h"
+#include "mlir/Analysis/Presburger/Fraction.h"
 #include "mlir/Analysis/Presburger/Utils.h"
 #include "llvm/Support/MathExtras.h"
 
 using namespace mlir;
 using namespace presburger;
 
-Matrix::Matrix(unsigned rows, unsigned columns, unsigned reservedRows,
+template <typename T> Matrix<T>::Matrix(unsigned rows, unsigned columns, unsigned reservedRows,
                unsigned reservedColumns)
     : nRows(rows), nColumns(columns),
       nReservedColumns(std::max(nColumns, reservedColumns)),

   data.reserve(std::max(nRows, reservedRows) * nReservedColumns);
 }
 
-Matrix Matrix::identity(unsigned dimension) {
+template <typename T> Matrix<T> Matrix<T>::identity(unsigned dimension) {
   Matrix matrix(dimension, dimension);
   for (unsigned i = 0; i < dimension; ++i)
     matrix(i, i) = 1;
   return matrix;
 }
 
-unsigned Matrix::getNumReservedRows() const {
+template <typename T> unsigned Matrix<T>::getNumReservedRows() const {
   return data.capacity() / nReservedColumns;
 }
 
-void Matrix::reserveRows(unsigned rows) {
+template <typename T> void Matrix<T>::reserveRows(unsigned rows) {
   data.reserve(rows * nReservedColumns);
 }
 
-unsigned Matrix::appendExtraRow() {
+template <typename T> unsigned Matrix<T>::appendExtraRow() {
   resizeVertically(nRows + 1);
   return nRows - 1;
 }
 
-unsigned Matrix::appendExtraRow(ArrayRef<MPInt> elems) {
+template <typename T> unsigned Matrix<T>::appendExtraRow(ArrayRef<T> elems) {
   assert(elems.size() == nColumns && "elems must match row length!");
   unsigned row = appendExtraRow();
   for (unsigned col = 0; col < nColumns; ++col)

   return row;
 }
 
-void Matrix::resizeHorizontally(unsigned newNColumns) {
+template <typename T> void Matrix<T>::resizeHorizontally(unsigned newNColumns) {
   if (newNColumns < nColumns)
     removeColumns(newNColumns, nColumns - newNColumns);
   if (newNColumns > nColumns)
     insertColumns(nColumns, newNColumns - nColumns);
 }
 
-void Matrix::resize(unsigned newNRows, unsigned newNColumns) {
+template <typename T> void Matrix<T>::resize(unsigned newNRows, unsigned newNColumns) {
   resizeHorizontally(newNColumns);
   resizeVertically(newNRows);
 }
 
-void Matrix::resizeVertically(unsigned newNRows) {
+template <typename T> void Matrix<T>::resizeVertically(unsigned newNRows) {
   nRows = newNRows;
   data.resize(nRows * nReservedColumns);
 }
 
-void Matrix::swapRows(unsigned row, unsigned otherRow) {
+template <typename T> void Matrix<T>::swapRows(unsigned row, unsigned otherRow) {
   assert((row < getNumRows() && otherRow < getNumRows()) &&
          "Given row out of bounds");
   if (row == otherRow)

     std::swap(at(row, col), at(otherRow, col));
 }
 
-void Matrix::swapColumns(unsigned column, unsigned otherColumn) {
+template <typename T> void Matrix<T>::swapColumns(unsigned column, unsigned otherColumn) {
   assert((column < getNumColumns() && otherColumn < getNumColumns()) &&
          "Given column out of bounds");
   if (column == otherColumn)

     std::swap(at(row, column), at(row, otherColumn));
 }
 
-MutableArrayRef<MPInt> Matrix::getRow(unsigned row) {
+template <typename T> MutableArrayRef<T> Matrix<T>::getRow(unsigned row) {
   return {&data[row * nReservedColumns], nColumns};
 }
 
-ArrayRef<MPInt> Matrix::getRow(unsigned row) const {
+template <typename T> ArrayRef<T> Matrix<T>::getRow(unsigned row) const {
   return {&data[row * nReservedColumns], nColumns};
 }
 
-void Matrix::setRow(unsigned row, ArrayRef<MPInt> elems) {
+template <typename T> void Matrix<T>::setRow(unsigned row, ArrayRef<T> elems) {
   assert(elems.size() == getNumColumns() &&
          "elems size must match row length!");
   for (unsigned i = 0, e = getNumColumns(); i < e; ++i)
     at(row, i) = elems[i];
 }
 
-void Matrix::insertColumn(unsigned pos) { insertColumns(pos, 1); }
-void Matrix::insertColumns(unsigned pos, unsigned count) {
+template <typename T> void Matrix<T>::insertColumn(unsigned pos) { insertColumns(pos, 1); }
+template <typename T> void Matrix<T>::insertColumns(unsigned pos, unsigned count) {
   if (count == 0)
     return;
   assert(pos <= nColumns);

     for (int ci = nReservedColumns - 1; ci >= 0; --ci) {
       unsigned r = ri;
       unsigned c = ci;
-      MPInt &dest = data[r * nReservedColumns + c];
+      T &dest = data[r * nReservedColumns + c];
       if (c >= nColumns) { // NOLINT
         // Out of bounds columns are zero-initialized. NOLINT because clang-tidy
         // complains about this branch being the same as the c >= pos one.

   }
 }
 
-void Matrix::removeColumn(unsigned pos) { removeColumns(pos, 1); }
-void Matrix::removeColumns(unsigned pos, unsigned count) {
+template <typename T> void Matrix<T>::removeColumn(unsigned pos) { removeColumns(pos, 1); }
+template <typename T> void Matrix<T>::removeColumns(unsigned pos, unsigned count) {
   if (count == 0)
     return;
   assert(pos + count - 1 < nColumns);

   nColumns -= count;
 }
 
-void Matrix::insertRow(unsigned pos) { insertRows(pos, 1); }
-void Matrix::insertRows(unsigned pos, unsigned count) {
+template <typename T> void Matrix<T>::insertRow(unsigned pos) { insertRows(pos, 1); }
+template <typename T> void Matrix<T>::insertRows(unsigned pos, unsigned count) {
   if (count == 0)
     return;
 

       at(r, c) = 0;
 }
 
-void Matrix::removeRow(unsigned pos) { removeRows(pos, 1); }
-void Matrix::removeRows(unsigned pos, unsigned count) {
+template <typename T> void Matrix<T>::removeRow(unsigned pos) { removeRows(pos, 1); }
+template <typename T> void Matrix<T>::removeRows(unsigned pos, unsigned count) {
   if (count == 0)
     return;
   assert(pos + count - 1 <= nRows);

   resizeVertically(nRows - count);
 }
 
-void Matrix::copyRow(unsigned sourceRow, unsigned targetRow) {
+template <typename T> void Matrix<T>::copyRow(unsigned sourceRow, unsigned targetRow) {
   if (sourceRow == targetRow)
     return;
   for (unsigned c = 0; c < nColumns; ++c)
     at(targetRow, c) = at(sourceRow, c);
 }
 
-void Matrix::fillRow(unsigned row, const MPInt &value) {
+template <typename T> void Matrix<T>::fillRow(unsigned row, const T &value) {
   for (unsigned col = 0; col < nColumns; ++col)
     at(row, col) = value;
 }
 
-void Matrix::addToRow(unsigned sourceRow, unsigned targetRow,
-                      const MPInt &scale) {
+template <typename T> void Matrix<T>::addToRow(unsigned sourceRow, unsigned targetRow,
+                      const T &scale) {
   addToRow(targetRow, getRow(sourceRow), scale);
 }
 
-void Matrix::addToRow(unsigned row, ArrayRef<MPInt> rowVec,
-                      const MPInt &scale) {
+template <typename T> void Matrix<T>::addToRow(unsigned row, ArrayRef<T> rowVec,
+                      const T &scale) {
   if (scale == 0)
     return;
   for (unsigned col = 0; col < nColumns; ++col)
-    at(row, col) += scale * rowVec[col];
+    at(row, col) = at(row, col) + scale * rowVec[col];
 }
 
-void Matrix::addToColumn(unsigned sourceColumn, unsigned targetColumn,
-                         const MPInt &scale) {
+template <typename T> void Matrix<T>::addToColumn(unsigned sourceColumn, unsigned targetColumn,
+                         const T &scale) {
   if (scale == 0)
     return;
   for (unsigned row = 0, e = getNumRows(); row < e; ++row)
-    at(row, targetColumn) += scale * at(row, sourceColumn);
+    at(row, targetColumn) = at(row, targetColumn) + scale * at(row, sourceColumn);
 }
 
-void Matrix::negateColumn(unsigned column) {
+template <typename T> void Matrix<T>::negateColumn(unsigned column) {
   for (unsigned row = 0, e = getNumRows(); row < e; ++row)
     at(row, column) = -at(row, column);
 }
 
-void Matrix::negateRow(unsigned row) {
+template <typename T> void Matrix<T>::negateRow(unsigned row) {
   for (unsigned column = 0, e = getNumColumns(); column < e; ++column)
     at(row, column) = -at(row, column);
 }
 
-MPInt Matrix::normalizeRow(unsigned row, unsigned cols) {
+template <> MPInt Matrix<MPInt>::normalizeRow(unsigned row, unsigned cols) {
   return normalizeRange(getRow(row).slice(0, cols));
 }
 
-MPInt Matrix::normalizeRow(unsigned row) {
+template <> MPInt Matrix<MPInt>::normalizeRow(unsigned row) {
   return normalizeRow(row, getNumColumns());
 }
 
-SmallVector<MPInt, 8> Matrix::preMultiplyWithRow(ArrayRef<MPInt> rowVec) const {
+template <typename T> SmallVector<T, 8> Matrix<T>::preMultiplyWithRow(ArrayRef<T> rowVec) const {
   assert(rowVec.size() == getNumRows() && "Invalid row vector dimension!");
 
-  SmallVector<MPInt, 8> result(getNumColumns(), MPInt(0));
+  SmallVector<T, 8> result(getNumColumns(), T(0));
   for (unsigned col = 0, e = getNumColumns(); col < e; ++col)
     for (unsigned i = 0, e = getNumRows(); i < e; ++i)
-      result[col] += rowVec[i] * at(i, col);
+      result[col] = result[col] + rowVec[i] * at(i, col);
   return result;
 }
 
-SmallVector<MPInt, 8>
-Matrix::postMultiplyWithColumn(ArrayRef<MPInt> colVec) const {
+template <typename T> SmallVector<T, 8>
+Matrix<T>::postMultiplyWithColumn(ArrayRef<T> colVec) const {
   assert(getNumColumns() == colVec.size() &&
          "Invalid column vector dimension!");
 
-  SmallVector<MPInt, 8> result(getNumRows(), MPInt(0));
+  SmallVector<T, 8> result(getNumRows(), T(0));
   for (unsigned row = 0, e = getNumRows(); row < e; row++)
     for (unsigned i = 0, e = getNumColumns(); i < e; i++)
-      result[row] += at(row, i) * colVec[i];
+      result[row] = result[row] + at(row, i) * colVec[i];
   return result;
 }
 

 /// 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) {
+static void modEntryColumnOperation(Matrix<MPInt> &m, unsigned row, unsigned sourceCol,
+                                    unsigned targetCol, Matrix<MPInt> &otherMatrix) {
   assert(m(row, sourceCol) != 0 && "Cannot divide by zero!");
   assert(m(row, sourceCol) > 0 && "Source must be positive!");
   MPInt ratio = -floorDiv(m(row, targetCol), m(row, sourceCol));

   otherMatrix.addToColumn(sourceCol, targetCol, ratio);
 }
 
-std::pair<Matrix, Matrix> Matrix::computeHermiteNormalForm() const {
+template <> std::pair<Matrix<MPInt>, Matrix<MPInt>> Matrix<MPInt>::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());
+  Matrix<MPInt> h = *this;
+  Matrix<MPInt> u = Matrix<MPInt>::identity(h.getNumColumns());
 
   unsigned echelonCol = 0;
   // Invariant: in all rows above row, all columns from echelonCol onwards

   return {h, u};
 }
 
-void Matrix::print(raw_ostream &os) const {
+template <typename T> void Matrix<T>::print(raw_ostream &os) const {
   for (unsigned row = 0; row < nRows; ++row) {
     for (unsigned column = 0; column < nColumns; ++column)
       os << at(row, column) << ' ';

   }
 }
 
-void Matrix::dump() const { print(llvm::errs()); }
+template <typename T> void Matrix<T>::dump() const { print(llvm::errs()); }
 
-bool Matrix::hasConsistentState() const {
+template <typename T> bool Matrix<T>::hasConsistentState() const {
   if (data.size() != nRows * nReservedColumns)
     return false;
   if (nColumns > nReservedColumns)

 #endif
   return true;
 }
+
+namespace mlir
+{
+  namespace presburger
+  {
+    template class Matrix<MPInt>;
+    template class Matrix<Fraction>;
+  }
+}
 No newline at end of file

mlir/lib/Analysis/Presburger/Simplex.cpp

 }
 
 void SymbolicLexSimplex::recordOutput(SymbolicLexOpt &result) const {
-  Matrix output(0, domainPoly.getNumVars() + 1);
+  Matrix<MPInt> output(0, domainPoly.getNumVars() + 1);
   output.reserveRows(result.lexopt.getNumOutputs());
   for (const Unknown &u : var) {
     if (u.isSymbol)

 ///
 /// When incrementing i, no cached f values get invalidated. However, the cached
 /// duals do get invalidated as the duals for the higher levels are different.
-void Simplex::reduceBasis(Matrix &basis, unsigned level) {
+void Simplex::reduceBasis(Matrix<MPInt> &basis, unsigned level) {
   const Fraction epsilon(3, 4);
 
   if (level == basis.getNumRows() - 1)

     return {};
 
   unsigned nDims = var.size();
-  Matrix basis = Matrix::identity(nDims);
+  Matrix<MPInt> basis = Matrix<MPInt>::identity(nDims);
 
   unsigned level = 0;
   // The snapshot just before constraining a direction to a value at each level.

mlir/unittests/Analysis/Presburger/LinearTransformTest.cpp

 using namespace mlir;
 using namespace presburger;
 
-void testColumnEchelonForm(const Matrix &m, unsigned expectedRank) {
+void testColumnEchelonForm(const Matrix<MPInt> &m, unsigned expectedRank) {
   unsigned lastAllowedNonZeroCol = 0;
   std::pair<unsigned, LinearTransform> result =
       LinearTransform::makeTransformToColumnEchelon(m);

 
 TEST(LinearTransformTest, transformToColumnEchelonTest) {
   // m1, m2, m3 are rank 1 matrices -- the first and second rows are identical.
-  Matrix m1(2, 2);
+  Matrix<MPInt> m1(2, 2);
   m1(0, 0) = 4;
   m1(0, 1) = -7;
   m1(1, 0) = 4;
   m1(1, 1) = -7;
   testColumnEchelonForm(m1, 1u);
 
-  Matrix m2(2, 2);
+  Matrix<MPInt> m2(2, 2);
   m2(0, 0) = -4;
   m2(0, 1) = 7;
   m2(1, 0) = 4;
   m2(1, 1) = -7;
   testColumnEchelonForm(m2, 1u);
 
-  Matrix m3(2, 2);
+  Matrix<MPInt> m3(2, 2);
   m3(0, 0) = -4;
   m3(0, 1) = -7;
   m3(1, 0) = -4;

   testColumnEchelonForm(m3, 1u);
 
   // m4, m5, m6 are rank 2 matrices -- the first and second rows are different.
-  Matrix m4(2, 2);
+  Matrix<MPInt> m4(2, 2);
   m4(0, 0) = 4;
   m4(0, 1) = -7;
   m4(1, 0) = -4;
   m4(1, 1) = -7;
   testColumnEchelonForm(m4, 2u);
 
-  Matrix m5(2, 2);
+  Matrix<MPInt> m5(2, 2);
   m5(0, 0) = -4;
   m5(0, 1) = 7;
   m5(1, 0) = 4;
   m5(1, 1) = 7;
   testColumnEchelonForm(m5, 2u);
 
-  Matrix m6(2, 2);
+  Matrix<MPInt> m6(2, 2);
   m6(0, 0) = -4;
   m6(0, 1) = -7;
   m6(1, 0) = 4;

mlir/unittests/Analysis/Presburger/MatrixTest.cpp

 //===----------------------------------------------------------------------===//
 
 #include "mlir/Analysis/Presburger/Matrix.h"
+#include "mlir/Analysis/Presburger/Fraction.h"
 #include "./Utils.h"
 #include <gmock/gmock.h>
 #include <gtest/gtest.h>

 using namespace presburger;
 
 TEST(MatrixTest, ReadWrite) {
-  Matrix mat(5, 5);
+  Matrix<MPInt> mat(5, 5);
   for (unsigned row = 0; row < 5; ++row)
     for (unsigned col = 0; col < 5; ++col)
       mat(row, col) = 10 * row + col;

 }
 
 TEST(MatrixTest, SwapColumns) {
-  Matrix mat(5, 5);
+  Matrix<MPInt> mat(5, 5);
   for (unsigned row = 0; row < 5; ++row)
     for (unsigned col = 0; col < 5; ++col)
       mat(row, col) = col == 3 ? 1 : 0;

 }
 
 TEST(MatrixTest, SwapRows) {
-  Matrix mat(5, 5);
+  Matrix<MPInt> mat(5, 5);
   for (unsigned row = 0; row < 5; ++row)
     for (unsigned col = 0; col < 5; ++col)
       mat(row, col) = row == 2 ? 1 : 0;

 }
 
 TEST(MatrixTest, resizeVertically) {
-  Matrix mat(5, 5);
+  Matrix<MPInt> mat(5, 5);
   EXPECT_EQ(mat.getNumRows(), 5u);
   EXPECT_EQ(mat.getNumColumns(), 5u);
   for (unsigned row = 0; row < 5; ++row)

 }
 
 TEST(MatrixTest, insertColumns) {
-  Matrix mat(5, 5, 5, 10);
+  Matrix<MPInt> mat(5, 5, 5, 10);
   EXPECT_EQ(mat.getNumRows(), 5u);
   EXPECT_EQ(mat.getNumColumns(), 5u);
   for (unsigned row = 0; row < 5; ++row)

 }
 
 TEST(MatrixTest, insertRows) {
-  Matrix mat(5, 5, 5, 10);
+  Matrix<MPInt> mat(5, 5, 5, 10);
   ASSERT_TRUE(mat.hasConsistentState());
   EXPECT_EQ(mat.getNumRows(), 5u);
   EXPECT_EQ(mat.getNumColumns(), 5u);

 }
 
 TEST(MatrixTest, resize) {
-  Matrix mat(5, 5);
+  Matrix<MPInt> mat(5, 5);
   EXPECT_EQ(mat.getNumRows(), 5u);
   EXPECT_EQ(mat.getNumColumns(), 5u);
   for (unsigned row = 0; row < 5; ++row)

       EXPECT_EQ(mat(row, col), row >= 3 || col >= 3 ? 0 : int(10 * row + col));
 }
 
-static void checkHermiteNormalForm(const Matrix &mat,
-                                   const Matrix &hermiteForm) {
+static void checkHermiteNormalForm(const Matrix<MPInt> &mat,
+                                   const Matrix<MPInt> &hermiteForm) {
   auto [h, u] = mat.computeHermiteNormalForm();
 
   for (unsigned row = 0; row < mat.getNumRows(); row++)

 
   {
     // 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}});
+    Matrix<MPInt> mat = makeIntMatrix(3, 3, {{2, 3, 6}, {3, 2, 3}, {17, 11, 16}});
+    Matrix<MPInt> hermiteForm = makeIntMatrix(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(
+    Matrix<MPInt> mat = makeIntMatrix(
         4, 4,
         {{-6, -1, -19, -20}, {0, 1, 0, 0}, {-5, 0, -15, -16}, {6, 0, 18, 19}});
-    Matrix hermiteForm = makeMatrix(
+    Matrix<MPInt> hermiteForm = makeIntMatrix(
         4, 4, {{1, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, 1, 0}, {0, 0, 0, 1}});
     checkHermiteNormalForm(mat, hermiteForm);
   }
 
   {
-    Matrix mat = makeMatrix(
+    Matrix<MPInt> mat = makeIntMatrix(
         4, 4, {{3, 3, 1, 4}, {0, 1, 0, 0}, {0, 0, 19, 16}, {0, 0, 0, 3}});
-    Matrix hermiteForm = makeMatrix(
+    Matrix<MPInt> hermiteForm = makeIntMatrix(
         4, 4, {{1, 0, 0, 0}, {0, 1, 0, 0}, {1, 0, 3, 0}, {18, 0, 54, 57}});
     checkHermiteNormalForm(mat, hermiteForm);
   }
 
   {
-    Matrix mat = makeMatrix(
+    Matrix<MPInt> mat = makeIntMatrix(
         4, 4, {{3, 3, 1, 4}, {0, 1, 0, 0}, {0, 0, 19, 16}, {0, 0, 0, 3}});
-    Matrix hermiteForm = makeMatrix(
+    Matrix<MPInt> hermiteForm = makeIntMatrix(
         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}});
+    Matrix<MPInt> mat =
+        makeIntMatrix(3, 5, {{0, 2, 0, 7, 1}, {-1, 0, 0, -3, 0}, {0, 4, 1, 0, 8}});
+    Matrix<MPInt> hermiteForm =
+        makeIntMatrix(3, 5, {{1, 0, 0, 0, 0}, {0, 1, 0, 0, 0}, {0, 0, 1, 0, 0}});
     checkHermiteNormalForm(mat, hermiteForm);
   }
 }

mlir/unittests/Analysis/Presburger/Parser.h

 
   // TODO: Add default constructor for MultiAffineFunction.
   MultiAffineFunction multiAff(PresburgerSpace::getRelationSpace(),
-                               Matrix(0, 1));
+                               Matrix<MPInt>(0, 1));
   if (getMultiAffineFunctionFromMap(parseAffineMap(str, &context), multiAff)
           .failed())
     llvm_unreachable(

mlir/unittests/Analysis/Presburger/Utils.h

 #include "mlir/Analysis/Presburger/PWMAFunction.h"
 #include "mlir/Analysis/Presburger/PresburgerRelation.h"
 #include "mlir/Analysis/Presburger/Simplex.h"
+#include "mlir/Analysis/Presburger/Matrix.h"
 #include "mlir/IR/MLIRContext.h"
 #include "mlir/Support/LLVM.h"
 

 namespace mlir {
 namespace presburger {
 
-inline Matrix makeMatrix(unsigned numRow, unsigned numColumns,
-                         ArrayRef<SmallVector<int64_t, 8>> matrix) {
-  Matrix results(numRow, numColumns);
+inline Matrix<MPInt> makeIntMatrix(unsigned numRow, unsigned numColumns,
+                         ArrayRef<SmallVector<int, 8>> matrix) {
+  Matrix<MPInt> results(numRow, numColumns);
+  assert(matrix.size() == numRow);
+  for (unsigned i = 0; i < numRow; ++i) {
+    assert(matrix[i].size() == numColumns &&
+           "Output expression has incorrect dimensionality!");
+    for (unsigned j = 0; j < numColumns; ++j)
+      results(i, j) = MPInt(matrix[i][j]);
+  }
+  return results;
+}
+
+inline Matrix<Fraction> makeFracMatrix(unsigned numRow, unsigned numColumns,
+                         ArrayRef<SmallVector<Fraction, 8>> matrix) {
+  Matrix<Fraction> results(numRow, numColumns);
   assert(matrix.size() == numRow);
   for (unsigned i = 0; i < numRow; ++i) {
     assert(matrix[i].size() == numColumns &&