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/IntegerPolyhedron.h"
11 
12 namespace mlir {
13 
14 LinearTransform::LinearTransform(Matrix &&oMatrix) : matrix(oMatrix) {}
15 LinearTransform::LinearTransform(const Matrix &oMatrix) : matrix(oMatrix) {}
16 
17 // Set M(row, targetCol) to its remainder on division by M(row, sourceCol)
18 // by subtracting from column targetCol an appropriate integer multiple of
19 // sourceCol. This brings M(row, targetCol) to the range [0, M(row, sourceCol)).
20 // Apply the same column operation to otherMatrix, with the same integer
21 // multiple.
22 static void modEntryColumnOperation(Matrix &m, unsigned row, unsigned sourceCol,
23                                     unsigned targetCol, Matrix &otherMatrix) {
24   assert(m(row, sourceCol) != 0 && "Cannot divide by zero!");
25   assert((m(row, sourceCol) > 0 && m(row, targetCol) > 0) &&
26          "Operands must be positive!");
27   int64_t ratio = m(row, targetCol) / m(row, sourceCol);
28   m.addToColumn(sourceCol, targetCol, -ratio);
29   otherMatrix.addToColumn(sourceCol, targetCol, -ratio);
30 }
31 
32 std::pair<unsigned, LinearTransform>
33 LinearTransform::makeTransformToColumnEchelon(Matrix m) {
34   // We start with an identity result matrix and perform operations on m
35   // until m is in column echelon form. We apply the same sequence of operations
36   // on resultMatrix to obtain a transform that takes m to column echelon
37   // form.
38   Matrix resultMatrix = Matrix::identity(m.getNumColumns());
39 
40   unsigned echelonCol = 0;
41   // Invariant: in all rows above row, all columns from echelonCol onwards
42   // are all zero elements. In an iteration, if the curent row has any non-zero
43   // elements echelonCol onwards, we bring one to echelonCol and use it to
44   // make all elements echelonCol + 1 onwards zero.
45   for (unsigned row = 0; row < m.getNumRows(); ++row) {
46     // Search row for a non-empty entry, starting at echelonCol.
47     unsigned nonZeroCol = echelonCol;
48     for (unsigned e = m.getNumColumns(); nonZeroCol < e; ++nonZeroCol) {
49       if (m(row, nonZeroCol) == 0)
50         continue;
51       break;
52     }
53 
54     // Continue to the next row with the same echelonCol if this row is all
55     // zeros from echelonCol onwards.
56     if (nonZeroCol == m.getNumColumns())
57       continue;
58 
59     // Bring the non-zero column to echelonCol. This doesn't affect rows
60     // above since they are all zero at these columns.
61     if (nonZeroCol != echelonCol) {
62       m.swapColumns(nonZeroCol, echelonCol);
63       resultMatrix.swapColumns(nonZeroCol, echelonCol);
64     }
65 
66     // Make m(row, echelonCol) non-negative.
67     if (m(row, echelonCol) < 0) {
68       m.negateColumn(echelonCol);
69       resultMatrix.negateColumn(echelonCol);
70     }
71 
72     // Make all the entries in row after echelonCol zero.
73     for (unsigned i = echelonCol + 1, e = m.getNumColumns(); i < e; ++i) {
74       // We make m(row, i) non-negative, and then apply the Euclidean GCD
75       // algorithm to (row, i) and (row, echelonCol). At the end, one of them
76       // has value equal to the gcd of the two entries, and the other is zero.
77 
78       if (m(row, i) < 0) {
79         m.negateColumn(i);
80         resultMatrix.negateColumn(i);
81       }
82 
83       unsigned targetCol = i, sourceCol = echelonCol;
84       // At every step, we set m(row, targetCol) %= m(row, sourceCol), and
85       // swap the indices sourceCol and targetCol. (not the columns themselves)
86       // This modulo is implemented as a subtraction
87       // m(row, targetCol) -= quotient * m(row, sourceCol),
88       // where quotient = floor(m(row, targetCol) / m(row, sourceCol)),
89       // which brings m(row, targetCol) to the range [0, m(row, sourceCol)).
90       //
91       // We are only allowed column operations; we perform the above
92       // for every row, i.e., the above subtraction is done as a column
93       // operation. This does not affect any rows above us since they are
94       // guaranteed to be zero at these columns.
95       while (m(row, targetCol) != 0 && m(row, sourceCol) != 0) {
96         modEntryColumnOperation(m, row, sourceCol, targetCol, resultMatrix);
97         std::swap(targetCol, sourceCol);
98       }
99 
100       // One of (row, echelonCol) and (row, i) is zero and the other is the gcd.
101       // Make it so that (row, echelonCol) holds the non-zero value.
102       if (m(row, echelonCol) == 0) {
103         m.swapColumns(i, echelonCol);
104         resultMatrix.swapColumns(i, echelonCol);
105       }
106     }
107 
108     ++echelonCol;
109   }
110 
111   return {echelonCol, LinearTransform(std::move(resultMatrix))};
112 }
113 
114 IntegerPolyhedron
115 LinearTransform::applyTo(const IntegerPolyhedron &poly) const {
116   IntegerPolyhedron result(poly.getNumIds());
117 
118   for (unsigned i = 0, e = poly.getNumEqualities(); i < e; ++i) {
119     ArrayRef<int64_t> eq = poly.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 = poly.getNumInequalities(); i < e; ++i) {
129     ArrayRef<int64_t> ineq = poly.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 
141 } // namespace mlir
142