MLIR  17.0.0git
Matrix.cpp
Go to the documentation of this file.
1 //===- Matrix.cpp - MLIR Matrix 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 
11 #include "llvm/Support/MathExtras.h"
12 
13 using namespace mlir;
14 using namespace presburger;
15 
16 Matrix::Matrix(unsigned rows, unsigned columns, unsigned reservedRows,
17  unsigned reservedColumns)
18  : nRows(rows), nColumns(columns),
19  nReservedColumns(std::max(nColumns, reservedColumns)),
20  data(nRows * nReservedColumns) {
21  data.reserve(std::max(nRows, reservedRows) * nReservedColumns);
22 }
23 
24 Matrix Matrix::identity(unsigned dimension) {
25  Matrix matrix(dimension, dimension);
26  for (unsigned i = 0; i < dimension; ++i)
27  matrix(i, i) = 1;
28  return matrix;
29 }
30 
31 unsigned Matrix::getNumReservedRows() const {
32  return data.capacity() / nReservedColumns;
33 }
34 
35 void Matrix::reserveRows(unsigned rows) {
36  data.reserve(rows * nReservedColumns);
37 }
38 
40  resizeVertically(nRows + 1);
41  return nRows - 1;
42 }
43 
45  assert(elems.size() == nColumns && "elems must match row length!");
46  unsigned row = appendExtraRow();
47  for (unsigned col = 0; col < nColumns; ++col)
48  at(row, col) = elems[col];
49  return row;
50 }
51 
52 void Matrix::resizeHorizontally(unsigned newNColumns) {
53  if (newNColumns < nColumns)
54  removeColumns(newNColumns, nColumns - newNColumns);
55  if (newNColumns > nColumns)
56  insertColumns(nColumns, newNColumns - nColumns);
57 }
58 
59 void Matrix::resize(unsigned newNRows, unsigned newNColumns) {
60  resizeHorizontally(newNColumns);
61  resizeVertically(newNRows);
62 }
63 
64 void Matrix::resizeVertically(unsigned newNRows) {
65  nRows = newNRows;
66  data.resize(nRows * nReservedColumns);
67 }
68 
69 void Matrix::swapRows(unsigned row, unsigned otherRow) {
70  assert((row < getNumRows() && otherRow < getNumRows()) &&
71  "Given row out of bounds");
72  if (row == otherRow)
73  return;
74  for (unsigned col = 0; col < nColumns; col++)
75  std::swap(at(row, col), at(otherRow, col));
76 }
77 
78 void Matrix::swapColumns(unsigned column, unsigned otherColumn) {
79  assert((column < getNumColumns() && otherColumn < getNumColumns()) &&
80  "Given column out of bounds");
81  if (column == otherColumn)
82  return;
83  for (unsigned row = 0; row < nRows; row++)
84  std::swap(at(row, column), at(row, otherColumn));
85 }
86 
88  return {&data[row * nReservedColumns], nColumns};
89 }
90 
91 ArrayRef<MPInt> Matrix::getRow(unsigned row) const {
92  return {&data[row * nReservedColumns], nColumns};
93 }
94 
95 void Matrix::setRow(unsigned row, ArrayRef<MPInt> elems) {
96  assert(elems.size() == getNumColumns() &&
97  "elems size must match row length!");
98  for (unsigned i = 0, e = getNumColumns(); i < e; ++i)
99  at(row, i) = elems[i];
100 }
101 
102 void Matrix::insertColumn(unsigned pos) { insertColumns(pos, 1); }
103 void Matrix::insertColumns(unsigned pos, unsigned count) {
104  if (count == 0)
105  return;
106  assert(pos <= nColumns);
107  unsigned oldNReservedColumns = nReservedColumns;
108  if (nColumns + count > nReservedColumns) {
109  nReservedColumns = llvm::NextPowerOf2(nColumns + count);
110  data.resize(nRows * nReservedColumns);
111  }
112  nColumns += count;
113 
114  for (int ri = nRows - 1; ri >= 0; --ri) {
115  for (int ci = nReservedColumns - 1; ci >= 0; --ci) {
116  unsigned r = ri;
117  unsigned c = ci;
118  MPInt &dest = data[r * nReservedColumns + c];
119  if (c >= nColumns) { // NOLINT
120  // Out of bounds columns are zero-initialized. NOLINT because clang-tidy
121  // complains about this branch being the same as the c >= pos one.
122  //
123  // TODO: this case can be skipped if the number of reserved columns
124  // didn't change.
125  dest = 0;
126  } else if (c >= pos + count) {
127  // Shift the data occuring after the inserted columns.
128  dest = data[r * oldNReservedColumns + c - count];
129  } else if (c >= pos) {
130  // The inserted columns are also zero-initialized.
131  dest = 0;
132  } else {
133  // The columns before the inserted columns stay at the same (row, col)
134  // but this corresponds to a different location in the linearized array
135  // if the number of reserved columns changed.
136  if (nReservedColumns == oldNReservedColumns)
137  break;
138  dest = data[r * oldNReservedColumns + c];
139  }
140  }
141  }
142 }
143 
144 void Matrix::removeColumn(unsigned pos) { removeColumns(pos, 1); }
145 void Matrix::removeColumns(unsigned pos, unsigned count) {
146  if (count == 0)
147  return;
148  assert(pos + count - 1 < nColumns);
149  for (unsigned r = 0; r < nRows; ++r) {
150  for (unsigned c = pos; c < nColumns - count; ++c)
151  at(r, c) = at(r, c + count);
152  for (unsigned c = nColumns - count; c < nColumns; ++c)
153  at(r, c) = 0;
154  }
155  nColumns -= count;
156 }
157 
158 void Matrix::insertRow(unsigned pos) { insertRows(pos, 1); }
159 void Matrix::insertRows(unsigned pos, unsigned count) {
160  if (count == 0)
161  return;
162 
163  assert(pos <= nRows);
164  resizeVertically(nRows + count);
165  for (int r = nRows - 1; r >= int(pos + count); --r)
166  copyRow(r - count, r);
167  for (int r = pos + count - 1; r >= int(pos); --r)
168  for (unsigned c = 0; c < nColumns; ++c)
169  at(r, c) = 0;
170 }
171 
172 void Matrix::removeRow(unsigned pos) { removeRows(pos, 1); }
173 void Matrix::removeRows(unsigned pos, unsigned count) {
174  if (count == 0)
175  return;
176  assert(pos + count - 1 <= nRows);
177  for (unsigned r = pos; r + count < nRows; ++r)
178  copyRow(r + count, r);
179  resizeVertically(nRows - count);
180 }
181 
182 void Matrix::copyRow(unsigned sourceRow, unsigned targetRow) {
183  if (sourceRow == targetRow)
184  return;
185  for (unsigned c = 0; c < nColumns; ++c)
186  at(targetRow, c) = at(sourceRow, c);
187 }
188 
189 void Matrix::fillRow(unsigned row, const MPInt &value) {
190  for (unsigned col = 0; col < nColumns; ++col)
191  at(row, col) = value;
192 }
193 
194 void Matrix::addToRow(unsigned sourceRow, unsigned targetRow,
195  const MPInt &scale) {
196  addToRow(targetRow, getRow(sourceRow), scale);
197 }
198 
199 void Matrix::addToRow(unsigned row, ArrayRef<MPInt> rowVec,
200  const MPInt &scale) {
201  if (scale == 0)
202  return;
203  for (unsigned col = 0; col < nColumns; ++col)
204  at(row, col) += scale * rowVec[col];
205 }
206 
207 void Matrix::addToColumn(unsigned sourceColumn, unsigned targetColumn,
208  const MPInt &scale) {
209  if (scale == 0)
210  return;
211  for (unsigned row = 0, e = getNumRows(); row < e; ++row)
212  at(row, targetColumn) += scale * at(row, sourceColumn);
213 }
214 
215 void Matrix::negateColumn(unsigned column) {
216  for (unsigned row = 0, e = getNumRows(); row < e; ++row)
217  at(row, column) = -at(row, column);
218 }
219 
220 void Matrix::negateRow(unsigned row) {
221  for (unsigned column = 0, e = getNumColumns(); column < e; ++column)
222  at(row, column) = -at(row, column);
223 }
224 
225 MPInt Matrix::normalizeRow(unsigned row, unsigned cols) {
226  return normalizeRange(getRow(row).slice(0, cols));
227 }
228 
230  return normalizeRow(row, getNumColumns());
231 }
232 
234  assert(rowVec.size() == getNumRows() && "Invalid row vector dimension!");
235 
237  for (unsigned col = 0, e = getNumColumns(); col < e; ++col)
238  for (unsigned i = 0, e = getNumRows(); i < e; ++i)
239  result[col] += rowVec[i] * at(i, col);
240  return result;
241 }
242 
245  assert(getNumColumns() == colVec.size() &&
246  "Invalid column vector dimension!");
247 
248  SmallVector<MPInt, 8> result(getNumRows(), MPInt(0));
249  for (unsigned row = 0, e = getNumRows(); row < e; row++)
250  for (unsigned i = 0, e = getNumColumns(); i < e; i++)
251  result[row] += at(row, i) * colVec[i];
252  return result;
253 }
254 
255 /// Set M(row, targetCol) to its remainder on division by M(row, sourceCol)
256 /// by subtracting from column targetCol an appropriate integer multiple of
257 /// sourceCol. This brings M(row, targetCol) to the range [0, M(row,
258 /// sourceCol)). Apply the same column operation to otherMatrix, with the same
259 /// integer multiple.
260 static void modEntryColumnOperation(Matrix &m, unsigned row, unsigned sourceCol,
261  unsigned targetCol, Matrix &otherMatrix) {
262  assert(m(row, sourceCol) != 0 && "Cannot divide by zero!");
263  assert(m(row, sourceCol) > 0 && "Source must be positive!");
264  MPInt ratio = -floorDiv(m(row, targetCol), m(row, sourceCol));
265  m.addToColumn(sourceCol, targetCol, ratio);
266  otherMatrix.addToColumn(sourceCol, targetCol, ratio);
267 }
268 
269 std::pair<Matrix, Matrix> Matrix::computeHermiteNormalForm() const {
270  // We start with u as an identity matrix and perform operations on h until h
271  // is in hermite normal form. We apply the same sequence of operations on u to
272  // obtain a transform that takes h to hermite normal form.
273  Matrix h = *this;
275 
276  unsigned echelonCol = 0;
277  // Invariant: in all rows above row, all columns from echelonCol onwards
278  // are all zero elements. In an iteration, if the curent row has any non-zero
279  // elements echelonCol onwards, we bring one to echelonCol and use it to
280  // make all elements echelonCol + 1 onwards zero.
281  for (unsigned row = 0; row < h.getNumRows(); ++row) {
282  // Search row for a non-empty entry, starting at echelonCol.
283  unsigned nonZeroCol = echelonCol;
284  for (unsigned e = h.getNumColumns(); nonZeroCol < e; ++nonZeroCol) {
285  if (h(row, nonZeroCol) == 0)
286  continue;
287  break;
288  }
289 
290  // Continue to the next row with the same echelonCol if this row is all
291  // zeros from echelonCol onwards.
292  if (nonZeroCol == h.getNumColumns())
293  continue;
294 
295  // Bring the non-zero column to echelonCol. This doesn't affect rows
296  // above since they are all zero at these columns.
297  if (nonZeroCol != echelonCol) {
298  h.swapColumns(nonZeroCol, echelonCol);
299  u.swapColumns(nonZeroCol, echelonCol);
300  }
301 
302  // Make h(row, echelonCol) non-negative.
303  if (h(row, echelonCol) < 0) {
304  h.negateColumn(echelonCol);
305  u.negateColumn(echelonCol);
306  }
307 
308  // Make all the entries in row after echelonCol zero.
309  for (unsigned i = echelonCol + 1, e = h.getNumColumns(); i < e; ++i) {
310  // We make h(row, i) non-negative, and then apply the Euclidean GCD
311  // algorithm to (row, i) and (row, echelonCol). At the end, one of them
312  // has value equal to the gcd of the two entries, and the other is zero.
313 
314  if (h(row, i) < 0) {
315  h.negateColumn(i);
316  u.negateColumn(i);
317  }
318 
319  unsigned targetCol = i, sourceCol = echelonCol;
320  // At every step, we set h(row, targetCol) %= h(row, sourceCol), and
321  // swap the indices sourceCol and targetCol. (not the columns themselves)
322  // This modulo is implemented as a subtraction
323  // h(row, targetCol) -= quotient * h(row, sourceCol),
324  // where quotient = floor(h(row, targetCol) / h(row, sourceCol)),
325  // which brings h(row, targetCol) to the range [0, h(row, sourceCol)).
326  //
327  // We are only allowed column operations; we perform the above
328  // for every row, i.e., the above subtraction is done as a column
329  // operation. This does not affect any rows above us since they are
330  // guaranteed to be zero at these columns.
331  while (h(row, targetCol) != 0 && h(row, sourceCol) != 0) {
332  modEntryColumnOperation(h, row, sourceCol, targetCol, u);
333  std::swap(targetCol, sourceCol);
334  }
335 
336  // One of (row, echelonCol) and (row, i) is zero and the other is the gcd.
337  // Make it so that (row, echelonCol) holds the non-zero value.
338  if (h(row, echelonCol) == 0) {
339  h.swapColumns(i, echelonCol);
340  u.swapColumns(i, echelonCol);
341  }
342  }
343 
344  // Make all entries before echelonCol non-negative and strictly smaller
345  // than the pivot entry.
346  for (unsigned i = 0; i < echelonCol; ++i)
347  modEntryColumnOperation(h, row, echelonCol, i, u);
348 
349  ++echelonCol;
350  }
351 
352  return {h, u};
353 }
354 
355 void Matrix::print(raw_ostream &os) const {
356  for (unsigned row = 0; row < nRows; ++row) {
357  for (unsigned column = 0; column < nColumns; ++column)
358  os << at(row, column) << ' ';
359  os << '\n';
360  }
361 }
362 
363 void Matrix::dump() const { print(llvm::errs()); }
364 
366  if (data.size() != nRows * nReservedColumns)
367  return false;
368  if (nColumns > nReservedColumns)
369  return false;
370  for (unsigned r = 0; r < nRows; ++r)
371  for (unsigned c = nColumns; c < nReservedColumns; ++c)
372  if (data[r * nReservedColumns + c] != 0)
373  return false;
374  return true;
375 }
static void modEntryColumnOperation(Matrix &m, unsigned row, unsigned sourceCol, unsigned targetCol, Matrix &otherMatrix)
Set M(row, targetCol) to its remainder on division by M(row, sourceCol) by subtracting from column ta...
Definition: Matrix.cpp:260
static Value max(ImplicitLocOpBuilder &builder, Value value, Value bound)
This class provides support for multi-precision arithmetic.
Definition: MPInt.h:87
This is a class to represent a resizable matrix.
Definition: Matrix.h:35
bool hasConsistentState() const
Return whether the Matrix is in a consistent state with all its invariants satisfied.
Definition: Matrix.cpp:365
void insertRows(unsigned pos, unsigned count)
Insert rows having positions pos, pos + 1, ...
Definition: Matrix.cpp:159
void swapColumns(unsigned column, unsigned otherColumn)
Swap the given columns.
Definition: Matrix.cpp:78
void removeColumn(unsigned pos)
Definition: Matrix.cpp:144
unsigned appendExtraRow()
Add an extra row at the bottom of the matrix and return its position.
Definition: Matrix.cpp:39
void print(raw_ostream &os) const
Print the matrix.
Definition: Matrix.cpp:355
void copyRow(unsigned sourceRow, unsigned targetRow)
Definition: Matrix.cpp:182
void insertColumn(unsigned pos)
Definition: Matrix.cpp:102
SmallVector< MPInt, 8 > postMultiplyWithColumn(ArrayRef< MPInt > colVec) const
The given vector is interpreted as a column vector v.
Definition: Matrix.cpp:244
static Matrix identity(unsigned dimension)
Return the identity matrix of the specified dimension.
Definition: Matrix.cpp:24
void removeColumns(unsigned pos, unsigned count)
Remove the columns having positions pos, pos + 1, ...
Definition: Matrix.cpp:145
void insertColumns(unsigned pos, unsigned count)
Insert columns having positions pos, pos + 1, ...
Definition: Matrix.cpp:103
MPInt normalizeRow(unsigned row, unsigned nCols)
Divide the first nCols of the specified row by their GCD.
Definition: Matrix.cpp:225
MPInt & at(unsigned row, unsigned column)
Access the element at the specified row and column.
Definition: Matrix.h:52
SmallVector< MPInt, 8 > preMultiplyWithRow(ArrayRef< MPInt > rowVec) const
The given vector is interpreted as a row vector v.
Definition: Matrix.cpp:233
void removeRow(unsigned pos)
Definition: Matrix.cpp:172
void resizeVertically(unsigned newNRows)
Definition: Matrix.cpp:64
void setRow(unsigned row, ArrayRef< MPInt > elems)
Set the specified row to elems.
Definition: Matrix.cpp:95
unsigned getNumReservedRows() const
Return the maximum number of rows/columns that can be added without incurring a reallocation.
Definition: Matrix.cpp:31
unsigned getNumColumns() const
Definition: Matrix.h:78
void fillRow(unsigned row, const MPInt &value)
Definition: Matrix.cpp:189
unsigned getNumRows() const
Definition: Matrix.h:76
void addToColumn(unsigned sourceColumn, unsigned targetColumn, const MPInt &scale)
Add scale multiples of the source column to the target column.
Definition: Matrix.cpp:207
void insertRow(unsigned pos)
Definition: Matrix.cpp:158
void swapRows(unsigned row, unsigned otherRow)
Swap the given rows.
Definition: Matrix.cpp:69
void resizeHorizontally(unsigned newNColumns)
Definition: Matrix.cpp:52
void addToRow(unsigned sourceRow, unsigned targetRow, const MPInt &scale)
Add scale multiples of the source row to the target row.
Definition: Matrix.cpp:194
std::pair< Matrix, Matrix > computeHermiteNormalForm() const
Given the current matrix M, returns the matrices H, U such that H is the column hermite normal form o...
Definition: Matrix.cpp:269
void reserveRows(unsigned rows)
Reserve enough space to resize to the specified number of rows without reallocations.
Definition: Matrix.cpp:35
void negateColumn(unsigned column)
Negate the specified column.
Definition: Matrix.cpp:215
void resize(unsigned newNRows, unsigned newNColumns)
Resize the matrix to the specified dimensions.
Definition: Matrix.cpp:59
MutableArrayRef< MPInt > getRow(unsigned row)
Get a [Mutable]ArrayRef corresponding to the specified row.
Definition: Matrix.cpp:87
void negateRow(unsigned row)
Negate the specified row.
Definition: Matrix.cpp:220
void removeRows(unsigned pos, unsigned count)
Remove the rows having positions pos, pos + 1, ...
Definition: Matrix.cpp:173
MPInt normalizeRange(MutableArrayRef< MPInt > range)
Divide the range by its gcd and return the gcd.
Definition: Utils.cpp:347
LLVM_ATTRIBUTE_ALWAYS_INLINE MPInt floorDiv(const MPInt &lhs, const MPInt &rhs)
Definition: MPInt.h:382
Include the generated interface declarations.