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 SmallVector<int64_t, 8>
115 LinearTransform::preMultiplyWithRow(ArrayRef<int64_t> rowVec) const {
116   assert(rowVec.size() == matrix.getNumRows() &&
117          "row vector dimension should match transform output dimension");
118 
119   SmallVector<int64_t, 8> result(matrix.getNumColumns(), 0);
120   for (unsigned col = 0, e = matrix.getNumColumns(); col < e; ++col)
121     for (unsigned i = 0, e = matrix.getNumRows(); i < e; ++i)
122       result[col] += rowVec[i] * matrix(i, col);
123   return result;
124 }
125 
126 SmallVector<int64_t, 8>
127 LinearTransform::postMultiplyWithColumn(ArrayRef<int64_t> colVec) const {
128   assert(matrix.getNumColumns() == colVec.size() &&
129          "column vector dimension should match transform input dimension");
130 
131   SmallVector<int64_t, 8> result(matrix.getNumRows(), 0);
132   for (unsigned row = 0, e = matrix.getNumRows(); row < e; row++)
133     for (unsigned i = 0, e = matrix.getNumColumns(); i < e; i++)
134       result[row] += matrix(row, i) * colVec[i];
135   return result;
136 }
137 
138 IntegerPolyhedron
139 LinearTransform::applyTo(const IntegerPolyhedron &poly) const {
140   IntegerPolyhedron result(poly.getNumIds());
141 
142   for (unsigned i = 0, e = poly.getNumEqualities(); i < e; ++i) {
143     ArrayRef<int64_t> eq = poly.getEquality(i);
144 
145     int64_t c = eq.back();
146 
147     SmallVector<int64_t, 8> newEq = preMultiplyWithRow(eq.drop_back());
148     newEq.push_back(c);
149     result.addEquality(newEq);
150   }
151 
152   for (unsigned i = 0, e = poly.getNumInequalities(); i < e; ++i) {
153     ArrayRef<int64_t> ineq = poly.getInequality(i);
154 
155     int64_t c = ineq.back();
156 
157     SmallVector<int64_t, 8> newIneq = preMultiplyWithRow(ineq.drop_back());
158     newIneq.push_back(c);
159     result.addInequality(newIneq);
160   }
161 
162   return result;
163 }
164 
165 } // namespace mlir
166