1 //===- LinearTransform.cpp - MLIR LinearTransform Class -------------------===//
2 //
3 // Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
4 // See https://llvm.org/LICENSE.txt for license information.
5 // SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
6 //
7 //===----------------------------------------------------------------------===//
8 
9 #include "mlir/Analysis/Presburger/LinearTransform.h"
10 #include "mlir/Analysis/Presburger/IntegerRelation.h"
11 
12 using namespace mlir;
13 using namespace presburger;
14 
15 LinearTransform::LinearTransform(Matrix &&oMatrix) : matrix(oMatrix) {}
16 LinearTransform::LinearTransform(const Matrix &oMatrix) : matrix(oMatrix) {}
17 
18 // Set M(row, targetCol) to its remainder on division by M(row, sourceCol)
19 // by subtracting from column targetCol an appropriate integer multiple of
20 // sourceCol. This brings M(row, targetCol) to the range [0, M(row, sourceCol)).
21 // Apply the same column operation to otherMatrix, with the same integer
22 // multiple.
23 static void modEntryColumnOperation(Matrix &m, unsigned row, unsigned sourceCol,
24                                     unsigned targetCol, Matrix &otherMatrix) {
25   assert(m(row, sourceCol) != 0 && "Cannot divide by zero!");
26   assert((m(row, sourceCol) > 0 && m(row, targetCol) > 0) &&
27          "Operands must be positive!");
28   int64_t ratio = m(row, targetCol) / m(row, sourceCol);
29   m.addToColumn(sourceCol, targetCol, -ratio);
30   otherMatrix.addToColumn(sourceCol, targetCol, -ratio);
31 }
32 
33 std::pair<unsigned, LinearTransform>
34 LinearTransform::makeTransformToColumnEchelon(Matrix m) {
35   // We start with an identity result matrix and perform operations on m
36   // until m is in column echelon form. We apply the same sequence of operations
37   // on resultMatrix to obtain a transform that takes m to column echelon
38   // form.
39   Matrix resultMatrix = Matrix::identity(m.getNumColumns());
40 
41   unsigned echelonCol = 0;
42   // Invariant: in all rows above row, all columns from echelonCol onwards
43   // are all zero elements. In an iteration, if the curent row has any non-zero
44   // elements echelonCol onwards, we bring one to echelonCol and use it to
45   // make all elements echelonCol + 1 onwards zero.
46   for (unsigned row = 0; row < m.getNumRows(); ++row) {
47     // Search row for a non-empty entry, starting at echelonCol.
48     unsigned nonZeroCol = echelonCol;
49     for (unsigned e = m.getNumColumns(); nonZeroCol < e; ++nonZeroCol) {
50       if (m(row, nonZeroCol) == 0)
51         continue;
52       break;
53     }
54 
55     // Continue to the next row with the same echelonCol if this row is all
56     // zeros from echelonCol onwards.
57     if (nonZeroCol == m.getNumColumns())
58       continue;
59 
60     // Bring the non-zero column to echelonCol. This doesn't affect rows
61     // above since they are all zero at these columns.
62     if (nonZeroCol != echelonCol) {
63       m.swapColumns(nonZeroCol, echelonCol);
64       resultMatrix.swapColumns(nonZeroCol, echelonCol);
65     }
66 
67     // Make m(row, echelonCol) non-negative.
68     if (m(row, echelonCol) < 0) {
69       m.negateColumn(echelonCol);
70       resultMatrix.negateColumn(echelonCol);
71     }
72 
73     // Make all the entries in row after echelonCol zero.
74     for (unsigned i = echelonCol + 1, e = m.getNumColumns(); i < e; ++i) {
75       // We make m(row, i) non-negative, and then apply the Euclidean GCD
76       // algorithm to (row, i) and (row, echelonCol). At the end, one of them
77       // has value equal to the gcd of the two entries, and the other is zero.
78 
79       if (m(row, i) < 0) {
80         m.negateColumn(i);
81         resultMatrix.negateColumn(i);
82       }
83 
84       unsigned targetCol = i, sourceCol = echelonCol;
85       // At every step, we set m(row, targetCol) %= m(row, sourceCol), and
86       // swap the indices sourceCol and targetCol. (not the columns themselves)
87       // This modulo is implemented as a subtraction
88       // m(row, targetCol) -= quotient * m(row, sourceCol),
89       // where quotient = floor(m(row, targetCol) / m(row, sourceCol)),
90       // which brings m(row, targetCol) to the range [0, m(row, sourceCol)).
91       //
92       // We are only allowed column operations; we perform the above
93       // for every row, i.e., the above subtraction is done as a column
94       // operation. This does not affect any rows above us since they are
95       // guaranteed to be zero at these columns.
96       while (m(row, targetCol) != 0 && m(row, sourceCol) != 0) {
97         modEntryColumnOperation(m, row, sourceCol, targetCol, resultMatrix);
98         std::swap(targetCol, sourceCol);
99       }
100 
101       // One of (row, echelonCol) and (row, i) is zero and the other is the gcd.
102       // Make it so that (row, echelonCol) holds the non-zero value.
103       if (m(row, echelonCol) == 0) {
104         m.swapColumns(i, echelonCol);
105         resultMatrix.swapColumns(i, echelonCol);
106       }
107     }
108 
109     ++echelonCol;
110   }
111 
112   return {echelonCol, LinearTransform(std::move(resultMatrix))};
113 }
114 
115 IntegerRelation LinearTransform::applyTo(const IntegerRelation &rel) const {
116   IntegerRelation result(rel.getSpace());
117 
118   for (unsigned i = 0, e = rel.getNumEqualities(); i < e; ++i) {
119     ArrayRef<int64_t> eq = rel.getEquality(i);
120 
121     int64_t c = eq.back();
122 
123     SmallVector<int64_t, 8> newEq = preMultiplyWithRow(eq.drop_back());
124     newEq.push_back(c);
125     result.addEquality(newEq);
126   }
127 
128   for (unsigned i = 0, e = rel.getNumInequalities(); i < e; ++i) {
129     ArrayRef<int64_t> ineq = rel.getInequality(i);
130 
131     int64_t c = ineq.back();
132 
133     SmallVector<int64_t, 8> newIneq = preMultiplyWithRow(ineq.drop_back());
134     newIneq.push_back(c);
135     result.addInequality(newIneq);
136   }
137 
138   return result;
139 }
140