[Mlir-commits] [mlir] 84ab06b - [MLIR][Presburger] Add Gram-Schmidt (#70843)
llvmlistbot at llvm.org
llvmlistbot at llvm.org
Wed Dec 13 00:28:51 PST 2023
Author: Abhinav271828
Date: 2023-12-13T08:28:47Z
New Revision: 84ab06ba2f37252ab40f84c39f17addce63eaa88
URL: https://github.com/llvm/llvm-project/commit/84ab06ba2f37252ab40f84c39f17addce63eaa88
DIFF: https://github.com/llvm/llvm-project/commit/84ab06ba2f37252ab40f84c39f17addce63eaa88.diff
LOG: [MLIR][Presburger] Add Gram-Schmidt (#70843)
Implement Gram-Schmidt orthogonalisation for the FracMatrix class.
This requires dotProduct, which has been added as a util.
Added:
Modified:
mlir/include/mlir/Analysis/Presburger/Matrix.h
mlir/include/mlir/Analysis/Presburger/Utils.h
mlir/lib/Analysis/Presburger/Matrix.cpp
mlir/lib/Analysis/Presburger/Utils.cpp
mlir/unittests/Analysis/Presburger/MatrixTest.cpp
Removed:
################################################################################
diff --git a/mlir/include/mlir/Analysis/Presburger/Matrix.h b/mlir/include/mlir/Analysis/Presburger/Matrix.h
index 4d9f13832e069..89fad85c0c337 100644
--- a/mlir/include/mlir/Analysis/Presburger/Matrix.h
+++ b/mlir/include/mlir/Analysis/Presburger/Matrix.h
@@ -265,6 +265,11 @@ class FracMatrix : public Matrix<Fraction> {
// does not exist, which happens iff det = 0.
// Assert-fails if the matrix is not square.
Fraction determinant(FracMatrix *inverse = nullptr) const;
+
+ // Computes the Gram-Schmidt orthogonalisation
+ // of the rows of matrix (cubic time).
+ // The rows of the matrix must be linearly independent.
+ FracMatrix gramSchmidt() const;
};
} // namespace presburger
diff --git a/mlir/include/mlir/Analysis/Presburger/Utils.h b/mlir/include/mlir/Analysis/Presburger/Utils.h
index a451ae8bf5572..20af0bfcd62ba 100644
--- a/mlir/include/mlir/Analysis/Presburger/Utils.h
+++ b/mlir/include/mlir/Analysis/Presburger/Utils.h
@@ -276,6 +276,11 @@ SmallVector<MPInt, 8> getNegatedCoeffs(ArrayRef<MPInt> coeffs);
/// a_1 x_1 + ... + a_n x_ + c < 0, i.e., -a_1 x_1 - ... - a_n x_ - c - 1 >= 0,
/// since all the variables are constrained to be integers.
SmallVector<MPInt, 8> getComplementIneq(ArrayRef<MPInt> ineq);
+
+/// Compute the dot product of two vectors.
+/// The vectors must have the same sizes.
+Fraction dotProduct(ArrayRef<Fraction> a, ArrayRef<Fraction> b);
+
} // namespace presburger
} // namespace mlir
diff --git a/mlir/lib/Analysis/Presburger/Matrix.cpp b/mlir/lib/Analysis/Presburger/Matrix.cpp
index fe461a1f95819..1fcc6d072b44b 100644
--- a/mlir/lib/Analysis/Presburger/Matrix.cpp
+++ b/mlir/lib/Analysis/Presburger/Matrix.cpp
@@ -466,8 +466,8 @@ FracMatrix FracMatrix::identity(unsigned dimension) {
FracMatrix::FracMatrix(IntMatrix m)
: FracMatrix(m.getNumRows(), m.getNumColumns()) {
- for (unsigned i = 0; i < m.getNumRows(); i++)
- for (unsigned j = 0; j < m.getNumColumns(); j++)
+ for (unsigned i = 0, r = m.getNumRows(); i < r; i++)
+ for (unsigned j = 0, c = m.getNumColumns(); j < c; j++)
this->at(i, j) = m.at(i, j);
}
@@ -554,4 +554,26 @@ Fraction FracMatrix::determinant(FracMatrix *inverse) const {
determinant *= m.at(i, i);
return determinant;
+}
+
+FracMatrix FracMatrix::gramSchmidt() const {
+ // Create a copy of the argument to store
+ // the orthogonalised version.
+ FracMatrix orth(*this);
+
+ // For each vector (row) in the matrix, subtract its unit
+ // projection along each of the previous vectors.
+ // This ensures that it has no component in the direction
+ // of any of the previous vectors.
+ for (unsigned i = 1, e = getNumRows(); i < e; i++) {
+ for (unsigned j = 0; j < i; j++) {
+ Fraction jNormSquared = dotProduct(orth.getRow(j), orth.getRow(j));
+ assert(jNormSquared != 0 && "some row became zero! Inputs to this "
+ "function must be linearly independent.");
+ Fraction projectionScale =
+ dotProduct(orth.getRow(i), orth.getRow(j)) / jNormSquared;
+ orth.addToRow(j, i, -projectionScale);
+ }
+ }
+ return orth;
}
\ No newline at end of file
diff --git a/mlir/lib/Analysis/Presburger/Utils.cpp b/mlir/lib/Analysis/Presburger/Utils.cpp
index 612712a1d9078..5e267d2045bc1 100644
--- a/mlir/lib/Analysis/Presburger/Utils.cpp
+++ b/mlir/lib/Analysis/Presburger/Utils.cpp
@@ -529,3 +529,12 @@ SmallVector<int64_t, 8> presburger::getInt64Vec(ArrayRef<MPInt> range) {
std::transform(range.begin(), range.end(), result.begin(), int64FromMPInt);
return result;
}
+
+Fraction presburger::dotProduct(ArrayRef<Fraction> a, ArrayRef<Fraction> b) {
+ assert(a.size() == b.size() &&
+ "dot product is only valid for vectors of equal sizes!");
+ Fraction sum = 0;
+ for (unsigned i = 0, e = a.size(); i < e; i++)
+ sum += a[i] * b[i];
+ return sum;
+}
\ No newline at end of file
diff --git a/mlir/unittests/Analysis/Presburger/MatrixTest.cpp b/mlir/unittests/Analysis/Presburger/MatrixTest.cpp
index d05b05e004c5c..508d4fa369c14 100644
--- a/mlir/unittests/Analysis/Presburger/MatrixTest.cpp
+++ b/mlir/unittests/Analysis/Presburger/MatrixTest.cpp
@@ -310,3 +310,71 @@ TEST(MatrixTest, intInverse) {
EXPECT_EQ(det, 0);
}
+
+TEST(MatrixTest, gramSchmidt) {
+ FracMatrix mat =
+ makeFracMatrix(3, 5,
+ {{Fraction(3, 1), Fraction(4, 1), Fraction(5, 1),
+ Fraction(12, 1), Fraction(19, 1)},
+ {Fraction(4, 1), Fraction(5, 1), Fraction(6, 1),
+ Fraction(13, 1), Fraction(20, 1)},
+ {Fraction(7, 1), Fraction(8, 1), Fraction(9, 1),
+ Fraction(16, 1), Fraction(24, 1)}});
+
+ FracMatrix gramSchmidt = makeFracMatrix(
+ 3, 5,
+ {{Fraction(3, 1), Fraction(4, 1), Fraction(5, 1), Fraction(12, 1),
+ Fraction(19, 1)},
+ {Fraction(142, 185), Fraction(383, 555), Fraction(68, 111),
+ Fraction(13, 185), Fraction(-262, 555)},
+ {Fraction(53, 463), Fraction(27, 463), Fraction(1, 463),
+ Fraction(-181, 463), Fraction(100, 463)}});
+
+ FracMatrix gs = mat.gramSchmidt();
+
+ EXPECT_EQ_FRAC_MATRIX(gs, gramSchmidt);
+ for (unsigned i = 0; i < 3u; i++)
+ for (unsigned j = i + 1; j < 3u; j++)
+ EXPECT_EQ(dotProduct(gs.getRow(i), gs.getRow(j)), 0);
+
+ mat = makeFracMatrix(3, 3,
+ {{Fraction(20, 1), Fraction(17, 1), Fraction(10, 1)},
+ {Fraction(20, 1), Fraction(18, 1), Fraction(6, 1)},
+ {Fraction(15, 1), Fraction(14, 1), Fraction(10, 1)}});
+
+ gramSchmidt = makeFracMatrix(
+ 3, 3,
+ {{Fraction(20, 1), Fraction(17, 1), Fraction(10, 1)},
+ {Fraction(460, 789), Fraction(1180, 789), Fraction(-2926, 789)},
+ {Fraction(-2925, 3221), Fraction(3000, 3221), Fraction(750, 3221)}});
+
+ gs = mat.gramSchmidt();
+
+ EXPECT_EQ_FRAC_MATRIX(gs, gramSchmidt);
+ for (unsigned i = 0; i < 3u; i++)
+ for (unsigned j = i + 1; j < 3u; j++)
+ EXPECT_EQ(dotProduct(gs.getRow(i), gs.getRow(j)), 0);
+
+ mat = makeFracMatrix(
+ 4, 4,
+ {{Fraction(1, 26), Fraction(13, 12), Fraction(34, 13), Fraction(7, 10)},
+ {Fraction(40, 23), Fraction(34, 1), Fraction(11, 19), Fraction(15, 1)},
+ {Fraction(21, 22), Fraction(10, 9), Fraction(4, 11), Fraction(14, 11)},
+ {Fraction(35, 22), Fraction(1, 15), Fraction(5, 8), Fraction(30, 1)}});
+
+ gs = mat.gramSchmidt();
+
+ // The integers involved are too big to construct the actual matrix.
+ // but we can check that the result is linearly independent.
+ ASSERT_FALSE(mat.determinant(nullptr) == 0);
+
+ for (unsigned i = 0; i < 4u; i++)
+ for (unsigned j = i + 1; j < 4u; j++)
+ EXPECT_EQ(dotProduct(gs.getRow(i), gs.getRow(j)), 0);
+
+ mat = FracMatrix::identity(/*dim=*/10);
+
+ gs = mat.gramSchmidt();
+
+ EXPECT_EQ_FRAC_MATRIX(gs, FracMatrix::identity(10));
+}
\ No newline at end of file
More information about the Mlir-commits
mailing list