MLIR  19.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 
13 #include "mlir/Support/LLVM.h"
14 #include "llvm/Support/MathExtras.h"
15 #include "llvm/Support/raw_ostream.h"
16 #include <algorithm>
17 #include <cassert>
18 #include <utility>
19 
20 using namespace mlir;
21 using namespace presburger;
22 
23 template <typename T>
24 Matrix<T>::Matrix(unsigned rows, unsigned columns, unsigned reservedRows,
25  unsigned reservedColumns)
26  : nRows(rows), nColumns(columns),
27  nReservedColumns(std::max(nColumns, reservedColumns)),
28  data(nRows * nReservedColumns) {
29  data.reserve(std::max(nRows, reservedRows) * nReservedColumns);
30 }
31 
32 /// We cannot use the default implementation of operator== as it compares
33 /// fields like `reservedColumns` etc., which are not part of the data.
34 template <typename T>
35 bool Matrix<T>::operator==(const Matrix<T> &m) const {
36  if (nRows != m.getNumRows())
37  return false;
38  if (nColumns != m.getNumColumns())
39  return false;
40 
41  for (unsigned i = 0; i < nRows; i++)
42  if (getRow(i) != m.getRow(i))
43  return false;
44 
45  return true;
46 }
47 
48 template <typename T>
49 Matrix<T> Matrix<T>::identity(unsigned dimension) {
50  Matrix matrix(dimension, dimension);
51  for (unsigned i = 0; i < dimension; ++i)
52  matrix(i, i) = 1;
53  return matrix;
54 }
55 
56 template <typename T>
58  return data.capacity() / nReservedColumns;
59 }
60 
61 template <typename T>
62 void Matrix<T>::reserveRows(unsigned rows) {
63  data.reserve(rows * nReservedColumns);
64 }
65 
66 template <typename T>
68  resizeVertically(nRows + 1);
69  return nRows - 1;
70 }
71 
72 template <typename T>
74  assert(elems.size() == nColumns && "elems must match row length!");
75  unsigned row = appendExtraRow();
76  for (unsigned col = 0; col < nColumns; ++col)
77  at(row, col) = elems[col];
78  return row;
79 }
80 
81 template <typename T>
83  Matrix<T> transp(nColumns, nRows);
84  for (unsigned row = 0; row < nRows; ++row)
85  for (unsigned col = 0; col < nColumns; ++col)
86  transp(col, row) = at(row, col);
87 
88  return transp;
89 }
90 
91 template <typename T>
92 void Matrix<T>::resizeHorizontally(unsigned newNColumns) {
93  if (newNColumns < nColumns)
94  removeColumns(newNColumns, nColumns - newNColumns);
95  if (newNColumns > nColumns)
96  insertColumns(nColumns, newNColumns - nColumns);
97 }
98 
99 template <typename T>
100 void Matrix<T>::resize(unsigned newNRows, unsigned newNColumns) {
101  resizeHorizontally(newNColumns);
102  resizeVertically(newNRows);
103 }
104 
105 template <typename T>
106 void Matrix<T>::resizeVertically(unsigned newNRows) {
107  nRows = newNRows;
108  data.resize(nRows * nReservedColumns);
109 }
110 
111 template <typename T>
112 void Matrix<T>::swapRows(unsigned row, unsigned otherRow) {
113  assert((row < getNumRows() && otherRow < getNumRows()) &&
114  "Given row out of bounds");
115  if (row == otherRow)
116  return;
117  for (unsigned col = 0; col < nColumns; col++)
118  std::swap(at(row, col), at(otherRow, col));
119 }
120 
121 template <typename T>
122 void Matrix<T>::swapColumns(unsigned column, unsigned otherColumn) {
123  assert((column < getNumColumns() && otherColumn < getNumColumns()) &&
124  "Given column out of bounds");
125  if (column == otherColumn)
126  return;
127  for (unsigned row = 0; row < nRows; row++)
128  std::swap(at(row, column), at(row, otherColumn));
129 }
130 
131 template <typename T>
133  return {&data[row * nReservedColumns], nColumns};
134 }
135 
136 template <typename T>
137 ArrayRef<T> Matrix<T>::getRow(unsigned row) const {
138  return {&data[row * nReservedColumns], nColumns};
139 }
140 
141 template <typename T>
142 void Matrix<T>::setRow(unsigned row, ArrayRef<T> elems) {
143  assert(elems.size() == getNumColumns() &&
144  "elems size must match row length!");
145  for (unsigned i = 0, e = getNumColumns(); i < e; ++i)
146  at(row, i) = elems[i];
147 }
148 
149 template <typename T>
150 void Matrix<T>::insertColumn(unsigned pos) {
151  insertColumns(pos, 1);
152 }
153 template <typename T>
154 void Matrix<T>::insertColumns(unsigned pos, unsigned count) {
155  if (count == 0)
156  return;
157  assert(pos <= nColumns);
158  unsigned oldNReservedColumns = nReservedColumns;
159  if (nColumns + count > nReservedColumns) {
160  nReservedColumns = llvm::NextPowerOf2(nColumns + count);
161  data.resize(nRows * nReservedColumns);
162  }
163  nColumns += count;
164 
165  for (int ri = nRows - 1; ri >= 0; --ri) {
166  for (int ci = nReservedColumns - 1; ci >= 0; --ci) {
167  unsigned r = ri;
168  unsigned c = ci;
169  T &dest = data[r * nReservedColumns + c];
170  if (c >= nColumns) { // NOLINT
171  // Out of bounds columns are zero-initialized. NOLINT because clang-tidy
172  // complains about this branch being the same as the c >= pos one.
173  //
174  // TODO: this case can be skipped if the number of reserved columns
175  // didn't change.
176  dest = 0;
177  } else if (c >= pos + count) {
178  // Shift the data occuring after the inserted columns.
179  dest = data[r * oldNReservedColumns + c - count];
180  } else if (c >= pos) {
181  // The inserted columns are also zero-initialized.
182  dest = 0;
183  } else {
184  // The columns before the inserted columns stay at the same (row, col)
185  // but this corresponds to a different location in the linearized array
186  // if the number of reserved columns changed.
187  if (nReservedColumns == oldNReservedColumns)
188  break;
189  dest = data[r * oldNReservedColumns + c];
190  }
191  }
192  }
193 }
194 
195 template <typename T>
196 void Matrix<T>::removeColumn(unsigned pos) {
197  removeColumns(pos, 1);
198 }
199 template <typename T>
200 void Matrix<T>::removeColumns(unsigned pos, unsigned count) {
201  if (count == 0)
202  return;
203  assert(pos + count - 1 < nColumns);
204  for (unsigned r = 0; r < nRows; ++r) {
205  for (unsigned c = pos; c < nColumns - count; ++c)
206  at(r, c) = at(r, c + count);
207  for (unsigned c = nColumns - count; c < nColumns; ++c)
208  at(r, c) = 0;
209  }
210  nColumns -= count;
211 }
212 
213 template <typename T>
214 void Matrix<T>::insertRow(unsigned pos) {
215  insertRows(pos, 1);
216 }
217 template <typename T>
218 void Matrix<T>::insertRows(unsigned pos, unsigned count) {
219  if (count == 0)
220  return;
221 
222  assert(pos <= nRows);
223  resizeVertically(nRows + count);
224  for (int r = nRows - 1; r >= int(pos + count); --r)
225  copyRow(r - count, r);
226  for (int r = pos + count - 1; r >= int(pos); --r)
227  for (unsigned c = 0; c < nColumns; ++c)
228  at(r, c) = 0;
229 }
230 
231 template <typename T>
232 void Matrix<T>::removeRow(unsigned pos) {
233  removeRows(pos, 1);
234 }
235 template <typename T>
236 void Matrix<T>::removeRows(unsigned pos, unsigned count) {
237  if (count == 0)
238  return;
239  assert(pos + count - 1 <= nRows);
240  for (unsigned r = pos; r + count < nRows; ++r)
241  copyRow(r + count, r);
242  resizeVertically(nRows - count);
243 }
244 
245 template <typename T>
246 void Matrix<T>::copyRow(unsigned sourceRow, unsigned targetRow) {
247  if (sourceRow == targetRow)
248  return;
249  for (unsigned c = 0; c < nColumns; ++c)
250  at(targetRow, c) = at(sourceRow, c);
251 }
252 
253 template <typename T>
254 void Matrix<T>::fillRow(unsigned row, const T &value) {
255  for (unsigned col = 0; col < nColumns; ++col)
256  at(row, col) = value;
257 }
258 
259 // moveColumns is implemented by moving the columns adjacent to the source range
260 // to their final position. When moving right (i.e. dstPos > srcPos), the range
261 // of the adjacent columns is [srcPos + num, dstPos + num). When moving left
262 // (i.e. dstPos < srcPos) the range of the adjacent columns is [dstPos, srcPos).
263 // First, zeroed out columns are inserted in the final positions of the adjacent
264 // columns. Then, the adjacent columns are moved to their final positions by
265 // swapping them with the zeroed columns. Finally, the now zeroed adjacent
266 // columns are deleted.
267 template <typename T>
268 void Matrix<T>::moveColumns(unsigned srcPos, unsigned num, unsigned dstPos) {
269  if (num == 0)
270  return;
271 
272  int offset = dstPos - srcPos;
273  if (offset == 0)
274  return;
275 
276  assert(srcPos + num <= getNumColumns() &&
277  "move source range exceeds matrix columns");
278  assert(dstPos + num <= getNumColumns() &&
279  "move destination range exceeds matrix columns");
280 
281  unsigned insertCount = offset > 0 ? offset : -offset;
282  unsigned finalAdjStart = offset > 0 ? srcPos : srcPos + num;
283  unsigned curAdjStart = offset > 0 ? srcPos + num : dstPos;
284  // TODO: This can be done using std::rotate.
285  // Insert new zero columns in the positions where the adjacent columns are to
286  // be moved.
287  insertColumns(finalAdjStart, insertCount);
288  // Update curAdjStart if insertion of new columns invalidates it.
289  if (finalAdjStart < curAdjStart)
290  curAdjStart += insertCount;
291 
292  // Swap the adjacent columns with inserted zero columns.
293  for (unsigned i = 0; i < insertCount; ++i)
294  swapColumns(finalAdjStart + i, curAdjStart + i);
295 
296  // Delete the now redundant zero columns.
297  removeColumns(curAdjStart, insertCount);
298 }
299 
300 template <typename T>
301 void Matrix<T>::addToRow(unsigned sourceRow, unsigned targetRow,
302  const T &scale) {
303  addToRow(targetRow, getRow(sourceRow), scale);
304 }
305 
306 template <typename T>
307 void Matrix<T>::addToRow(unsigned row, ArrayRef<T> rowVec, const T &scale) {
308  if (scale == 0)
309  return;
310  for (unsigned col = 0; col < nColumns; ++col)
311  at(row, col) += scale * rowVec[col];
312 }
313 
314 template <typename T>
315 void Matrix<T>::scaleRow(unsigned row, const T &scale) {
316  for (unsigned col = 0; col < nColumns; ++col)
317  at(row, col) *= scale;
318 }
319 
320 template <typename T>
321 void Matrix<T>::addToColumn(unsigned sourceColumn, unsigned targetColumn,
322  const T &scale) {
323  if (scale == 0)
324  return;
325  for (unsigned row = 0, e = getNumRows(); row < e; ++row)
326  at(row, targetColumn) += scale * at(row, sourceColumn);
327 }
328 
329 template <typename T>
330 void Matrix<T>::negateColumn(unsigned column) {
331  for (unsigned row = 0, e = getNumRows(); row < e; ++row)
332  at(row, column) = -at(row, column);
333 }
334 
335 template <typename T>
336 void Matrix<T>::negateRow(unsigned row) {
337  for (unsigned column = 0, e = getNumColumns(); column < e; ++column)
338  at(row, column) = -at(row, column);
339 }
340 
341 template <typename T>
343  for (unsigned row = 0; row < nRows; ++row)
344  negateRow(row);
345 }
346 
347 template <typename T>
349  assert(rowVec.size() == getNumRows() && "Invalid row vector dimension!");
350 
351  SmallVector<T, 8> result(getNumColumns(), T(0));
352  for (unsigned col = 0, e = getNumColumns(); col < e; ++col)
353  for (unsigned i = 0, e = getNumRows(); i < e; ++i)
354  result[col] += rowVec[i] * at(i, col);
355  return result;
356 }
357 
358 template <typename T>
360  assert(getNumColumns() == colVec.size() &&
361  "Invalid column vector dimension!");
362 
363  SmallVector<T, 8> result(getNumRows(), T(0));
364  for (unsigned row = 0, e = getNumRows(); row < e; row++)
365  for (unsigned i = 0, e = getNumColumns(); i < e; i++)
366  result[row] += at(row, i) * colVec[i];
367  return result;
368 }
369 
370 /// Set M(row, targetCol) to its remainder on division by M(row, sourceCol)
371 /// by subtracting from column targetCol an appropriate integer multiple of
372 /// sourceCol. This brings M(row, targetCol) to the range [0, M(row,
373 /// sourceCol)). Apply the same column operation to otherMatrix, with the same
374 /// integer multiple.
375 static void modEntryColumnOperation(Matrix<MPInt> &m, unsigned row,
376  unsigned sourceCol, unsigned targetCol,
377  Matrix<MPInt> &otherMatrix) {
378  assert(m(row, sourceCol) != 0 && "Cannot divide by zero!");
379  assert(m(row, sourceCol) > 0 && "Source must be positive!");
380  MPInt ratio = -floorDiv(m(row, targetCol), m(row, sourceCol));
381  m.addToColumn(sourceCol, targetCol, ratio);
382  otherMatrix.addToColumn(sourceCol, targetCol, ratio);
383 }
384 
385 template <typename T>
386 Matrix<T> Matrix<T>::getSubMatrix(unsigned fromRow, unsigned toRow,
387  unsigned fromColumn,
388  unsigned toColumn) const {
389  assert(fromRow <= toRow && "end of row range must be after beginning!");
390  assert(toRow < nRows && "end of row range out of bounds!");
391  assert(fromColumn <= toColumn &&
392  "end of column range must be after beginning!");
393  assert(toColumn < nColumns && "end of column range out of bounds!");
394  Matrix<T> subMatrix(toRow - fromRow + 1, toColumn - fromColumn + 1);
395  for (unsigned i = fromRow; i <= toRow; ++i)
396  for (unsigned j = fromColumn; j <= toColumn; ++j)
397  subMatrix(i - fromRow, j - fromColumn) = at(i, j);
398  return subMatrix;
399 }
400 
401 template <typename T>
402 void Matrix<T>::print(raw_ostream &os) const {
403  for (unsigned row = 0; row < nRows; ++row) {
404  for (unsigned column = 0; column < nColumns; ++column)
405  os << at(row, column) << ' ';
406  os << '\n';
407  }
408 }
409 
410 /// We iterate over the `indicator` bitset, checking each bit. If a bit is 1,
411 /// we append it to one matrix, and if it is zero, we append it to the other.
412 template <typename T>
413 std::pair<Matrix<T>, Matrix<T>>
415  Matrix<T> rowsForOne(0, nColumns), rowsForZero(0, nColumns);
416  for (unsigned i = 0; i < nRows; i++) {
417  if (indicator[i] == 1)
418  rowsForOne.appendExtraRow(getRow(i));
419  else
420  rowsForZero.appendExtraRow(getRow(i));
421  }
422  return {rowsForOne, rowsForZero};
423 }
424 
425 template <typename T>
426 void Matrix<T>::dump() const {
427  print(llvm::errs());
428 }
429 
430 template <typename T>
432  if (data.size() != nRows * nReservedColumns)
433  return false;
434  if (nColumns > nReservedColumns)
435  return false;
436 #ifdef EXPENSIVE_CHECKS
437  for (unsigned r = 0; r < nRows; ++r)
438  for (unsigned c = nColumns; c < nReservedColumns; ++c)
439  if (data[r * nReservedColumns + c] != 0)
440  return false;
441 #endif
442  return true;
443 }
444 
445 namespace mlir {
446 namespace presburger {
447 template class Matrix<MPInt>;
448 template class Matrix<Fraction>;
449 } // namespace presburger
450 } // namespace mlir
451 
452 IntMatrix IntMatrix::identity(unsigned dimension) {
453  IntMatrix matrix(dimension, dimension);
454  for (unsigned i = 0; i < dimension; ++i)
455  matrix(i, i) = 1;
456  return matrix;
457 }
458 
459 std::pair<IntMatrix, IntMatrix> IntMatrix::computeHermiteNormalForm() const {
460  // We start with u as an identity matrix and perform operations on h until h
461  // is in hermite normal form. We apply the same sequence of operations on u to
462  // obtain a transform that takes h to hermite normal form.
463  IntMatrix h = *this;
465 
466  unsigned echelonCol = 0;
467  // Invariant: in all rows above row, all columns from echelonCol onwards
468  // are all zero elements. In an iteration, if the curent row has any non-zero
469  // elements echelonCol onwards, we bring one to echelonCol and use it to
470  // make all elements echelonCol + 1 onwards zero.
471  for (unsigned row = 0; row < h.getNumRows(); ++row) {
472  // Search row for a non-empty entry, starting at echelonCol.
473  unsigned nonZeroCol = echelonCol;
474  for (unsigned e = h.getNumColumns(); nonZeroCol < e; ++nonZeroCol) {
475  if (h(row, nonZeroCol) == 0)
476  continue;
477  break;
478  }
479 
480  // Continue to the next row with the same echelonCol if this row is all
481  // zeros from echelonCol onwards.
482  if (nonZeroCol == h.getNumColumns())
483  continue;
484 
485  // Bring the non-zero column to echelonCol. This doesn't affect rows
486  // above since they are all zero at these columns.
487  if (nonZeroCol != echelonCol) {
488  h.swapColumns(nonZeroCol, echelonCol);
489  u.swapColumns(nonZeroCol, echelonCol);
490  }
491 
492  // Make h(row, echelonCol) non-negative.
493  if (h(row, echelonCol) < 0) {
494  h.negateColumn(echelonCol);
495  u.negateColumn(echelonCol);
496  }
497 
498  // Make all the entries in row after echelonCol zero.
499  for (unsigned i = echelonCol + 1, e = h.getNumColumns(); i < e; ++i) {
500  // We make h(row, i) non-negative, and then apply the Euclidean GCD
501  // algorithm to (row, i) and (row, echelonCol). At the end, one of them
502  // has value equal to the gcd of the two entries, and the other is zero.
503 
504  if (h(row, i) < 0) {
505  h.negateColumn(i);
506  u.negateColumn(i);
507  }
508 
509  unsigned targetCol = i, sourceCol = echelonCol;
510  // At every step, we set h(row, targetCol) %= h(row, sourceCol), and
511  // swap the indices sourceCol and targetCol. (not the columns themselves)
512  // This modulo is implemented as a subtraction
513  // h(row, targetCol) -= quotient * h(row, sourceCol),
514  // where quotient = floor(h(row, targetCol) / h(row, sourceCol)),
515  // which brings h(row, targetCol) to the range [0, h(row, sourceCol)).
516  //
517  // We are only allowed column operations; we perform the above
518  // for every row, i.e., the above subtraction is done as a column
519  // operation. This does not affect any rows above us since they are
520  // guaranteed to be zero at these columns.
521  while (h(row, targetCol) != 0 && h(row, sourceCol) != 0) {
522  modEntryColumnOperation(h, row, sourceCol, targetCol, u);
523  std::swap(targetCol, sourceCol);
524  }
525 
526  // One of (row, echelonCol) and (row, i) is zero and the other is the gcd.
527  // Make it so that (row, echelonCol) holds the non-zero value.
528  if (h(row, echelonCol) == 0) {
529  h.swapColumns(i, echelonCol);
530  u.swapColumns(i, echelonCol);
531  }
532  }
533 
534  // Make all entries before echelonCol non-negative and strictly smaller
535  // than the pivot entry.
536  for (unsigned i = 0; i < echelonCol; ++i)
537  modEntryColumnOperation(h, row, echelonCol, i, u);
538 
539  ++echelonCol;
540  }
541 
542  return {h, u};
543 }
544 
545 MPInt IntMatrix::normalizeRow(unsigned row, unsigned cols) {
546  return normalizeRange(getRow(row).slice(0, cols));
547 }
548 
550  return normalizeRow(row, getNumColumns());
551 }
552 
554  assert(nRows == nColumns &&
555  "determinant can only be calculated for square matrices!");
556 
557  FracMatrix m(*this);
558 
559  FracMatrix fracInverse(nRows, nColumns);
560  MPInt detM = m.determinant(&fracInverse).getAsInteger();
561 
562  if (detM == 0)
563  return MPInt(0);
564 
565  if (!inverse)
566  return detM;
567 
568  *inverse = IntMatrix(nRows, nColumns);
569  for (unsigned i = 0; i < nRows; i++)
570  for (unsigned j = 0; j < nColumns; j++)
571  inverse->at(i, j) = (fracInverse.at(i, j) * detM).getAsInteger();
572 
573  return detM;
574 }
575 
576 FracMatrix FracMatrix::identity(unsigned dimension) {
577  return Matrix::identity(dimension);
578 }
579 
581  : FracMatrix(m.getNumRows(), m.getNumColumns()) {
582  for (unsigned i = 0, r = m.getNumRows(); i < r; i++)
583  for (unsigned j = 0, c = m.getNumColumns(); j < c; j++)
584  this->at(i, j) = m.at(i, j);
585 }
586 
588  assert(nRows == nColumns &&
589  "determinant can only be calculated for square matrices!");
590 
591  FracMatrix m(*this);
592  FracMatrix tempInv(nRows, nColumns);
593  if (inverse)
594  tempInv = FracMatrix::identity(nRows);
595 
596  Fraction a, b;
597  // Make the matrix into upper triangular form using
598  // gaussian elimination with row operations.
599  // If inverse is required, we apply more operations
600  // to turn the matrix into diagonal form. We apply
601  // the same operations to the inverse matrix,
602  // which is initially identity.
603  // Either way, the product of the diagonal elements
604  // is then the determinant.
605  for (unsigned i = 0; i < nRows; i++) {
606  if (m(i, i) == 0)
607  // First ensure that the diagonal
608  // element is nonzero, by swapping
609  // it with a nonzero row.
610  for (unsigned j = i + 1; j < nRows; j++) {
611  if (m(j, i) != 0) {
612  m.swapRows(j, i);
613  if (inverse)
614  tempInv.swapRows(j, i);
615  break;
616  }
617  }
618 
619  b = m.at(i, i);
620  if (b == 0)
621  return 0;
622 
623  // Set all elements above the
624  // diagonal to zero.
625  if (inverse) {
626  for (unsigned j = 0; j < i; j++) {
627  if (m.at(j, i) == 0)
628  continue;
629  a = m.at(j, i);
630  // Set element (j, i) to zero
631  // by subtracting the ith row,
632  // appropriately scaled.
633  m.addToRow(i, j, -a / b);
634  tempInv.addToRow(i, j, -a / b);
635  }
636  }
637 
638  // Set all elements below the
639  // diagonal to zero.
640  for (unsigned j = i + 1; j < nRows; j++) {
641  if (m.at(j, i) == 0)
642  continue;
643  a = m.at(j, i);
644  // Set element (j, i) to zero
645  // by subtracting the ith row,
646  // appropriately scaled.
647  m.addToRow(i, j, -a / b);
648  if (inverse)
649  tempInv.addToRow(i, j, -a / b);
650  }
651  }
652 
653  // Now only diagonal elements of m are nonzero, but they are
654  // not necessarily 1. To get the true inverse, we should
655  // normalize them and apply the same scale to the inverse matrix.
656  // For efficiency we skip scaling m and just scale tempInv appropriately.
657  if (inverse) {
658  for (unsigned i = 0; i < nRows; i++)
659  for (unsigned j = 0; j < nRows; j++)
660  tempInv.at(i, j) = tempInv.at(i, j) / m(i, i);
661 
662  *inverse = std::move(tempInv);
663  }
664 
665  Fraction determinant = 1;
666  for (unsigned i = 0; i < nRows; i++)
667  determinant *= m.at(i, i);
668 
669  return determinant;
670 }
671 
673  // Create a copy of the argument to store
674  // the orthogonalised version.
675  FracMatrix orth(*this);
676 
677  // For each vector (row) in the matrix, subtract its unit
678  // projection along each of the previous vectors.
679  // This ensures that it has no component in the direction
680  // of any of the previous vectors.
681  for (unsigned i = 1, e = getNumRows(); i < e; i++) {
682  for (unsigned j = 0; j < i; j++) {
683  Fraction jNormSquared = dotProduct(orth.getRow(j), orth.getRow(j));
684  assert(jNormSquared != 0 && "some row became zero! Inputs to this "
685  "function must be linearly independent.");
686  Fraction projectionScale =
687  dotProduct(orth.getRow(i), orth.getRow(j)) / jNormSquared;
688  orth.addToRow(j, i, -projectionScale);
689  }
690  }
691  return orth;
692 }
693 
694 // Convert the matrix, interpreted (row-wise) as a basis
695 // to an LLL-reduced basis.
696 //
697 // This is an implementation of the algorithm described in
698 // "Factoring polynomials with rational coefficients" by
699 // A. K. Lenstra, H. W. Lenstra Jr., L. Lovasz.
700 //
701 // Let {b_1, ..., b_n} be the current basis and
702 // {b_1*, ..., b_n*} be the Gram-Schmidt orthogonalised
703 // basis (unnormalized).
704 // Define the Gram-Schmidt coefficients μ_ij as
705 // (b_i • b_j*) / (b_j* • b_j*), where (•) represents the inner product.
706 //
707 // We iterate starting from the second row to the last row.
708 //
709 // For the kth row, we first check μ_kj for all rows j < k.
710 // We subtract b_j (scaled by the integer nearest to μ_kj)
711 // from b_k.
712 //
713 // Now, we update k.
714 // If b_k and b_{k-1} satisfy the Lovasz condition
715 // |b_k|^2 ≥ (δ - μ_k{k-1}^2) |b_{k-1}|^2,
716 // we are done and we increment k.
717 // Otherwise, we swap b_k and b_{k-1} and decrement k.
718 //
719 // We repeat this until k = n and return.
721  MPInt nearest;
722  Fraction mu;
723 
724  // `gsOrth` holds the Gram-Schmidt orthogonalisation
725  // of the matrix at all times. It is recomputed every
726  // time the matrix is modified during the algorithm.
727  // This is naive and can be optimised.
728  FracMatrix gsOrth = gramSchmidt();
729 
730  // We start from the second row.
731  unsigned k = 1;
732  while (k < getNumRows()) {
733  for (unsigned j = k - 1; j < k; j--) {
734  // Compute the Gram-Schmidt coefficient μ_jk.
735  mu = dotProduct(getRow(k), gsOrth.getRow(j)) /
736  dotProduct(gsOrth.getRow(j), gsOrth.getRow(j));
737  nearest = round(mu);
738  // Subtract b_j scaled by the integer nearest to μ_jk from b_k.
739  addToRow(k, getRow(j), -Fraction(nearest, 1));
740  gsOrth = gramSchmidt(); // Update orthogonalization.
741  }
742  mu = dotProduct(getRow(k), gsOrth.getRow(k - 1)) /
743  dotProduct(gsOrth.getRow(k - 1), gsOrth.getRow(k - 1));
744  // Check the Lovasz condition for b_k and b_{k-1}.
745  if (dotProduct(gsOrth.getRow(k), gsOrth.getRow(k)) >
746  (delta - mu * mu) *
747  dotProduct(gsOrth.getRow(k - 1), gsOrth.getRow(k - 1))) {
748  // If it is satisfied, proceed to the next k.
749  k += 1;
750  } else {
751  // If it is not satisfied, decrement k (without
752  // going beyond the second row).
753  swapRows(k, k - 1);
754  gsOrth = gramSchmidt(); // Update orthogonalization.
755  k = k > 1 ? k - 1 : 1;
756  }
757  }
758 }
759 
761  unsigned numRows = getNumRows();
762  unsigned numColumns = getNumColumns();
763  IntMatrix normalized(numRows, numColumns);
764 
765  MPInt lcmDenoms = MPInt(1);
766  for (unsigned i = 0; i < numRows; i++) {
767  // For a row, first compute the LCM of the denominators.
768  for (unsigned j = 0; j < numColumns; j++)
769  lcmDenoms = lcm(lcmDenoms, at(i, j).den);
770  // Then, multiply by it throughout and convert to integers.
771  for (unsigned j = 0; j < numColumns; j++)
772  normalized(i, j) = (at(i, j) * lcmDenoms).getAsInteger();
773  }
774  return normalized;
775 }
static void modEntryColumnOperation(Matrix< MPInt > &m, unsigned row, unsigned sourceCol, unsigned targetCol, Matrix< MPInt > &otherMatrix)
Set M(row, targetCol) to its remainder on division by M(row, sourceCol) by subtracting from column ta...
Definition: Matrix.cpp:375
static Value max(ImplicitLocOpBuilder &builder, Value value, Value bound)
static void print(spirv::VerCapExtAttr triple, DialectAsmPrinter &printer)
FracMatrix gramSchmidt() const
Definition: Matrix.cpp:672
IntMatrix normalizeRows() const
Definition: Matrix.cpp:760
static FracMatrix identity(unsigned dimension)
Return the identity matrix of the specified dimension.
Definition: Matrix.cpp:576
FracMatrix(unsigned rows, unsigned columns, unsigned reservedRows=0, unsigned reservedColumns=0)
Definition: Matrix.h:294
void LLL(Fraction delta)
Definition: Matrix.cpp:720
Fraction determinant(FracMatrix *inverse=nullptr) const
Definition: Matrix.cpp:587
MPInt normalizeRow(unsigned row, unsigned nCols)
Divide the first nCols of the specified row by their GCD.
Definition: Matrix.cpp:545
static IntMatrix identity(unsigned dimension)
Return the identity matrix of the specified dimension.
Definition: Matrix.cpp:452
std::pair< IntMatrix, IntMatrix > 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:459
MPInt determinant(IntMatrix *inverse=nullptr) const
Definition: Matrix.cpp:553
IntMatrix(unsigned rows, unsigned columns, unsigned reservedRows=0, unsigned reservedColumns=0)
Definition: Matrix.h:252
This class provides support for multi-precision arithmetic.
Definition: MPInt.h:87
This is a class to represent a resizable matrix.
Definition: Matrix.h:41
void moveColumns(unsigned srcPos, unsigned num, unsigned dstPos)
Move the columns in the source range [srcPos, srcPos + num) to the specified destination [dstPos,...
Definition: Matrix.cpp:268
bool hasConsistentState() const
Return whether the Matrix is in a consistent state with all its invariants satisfied.
Definition: Matrix.cpp:431
void insertRows(unsigned pos, unsigned count)
Insert rows having positions pos, pos + 1, ...
Definition: Matrix.cpp:218
unsigned getNumRows() const
Definition: Matrix.h:85
void swapColumns(unsigned column, unsigned otherColumn)
Swap the given columns.
Definition: Matrix.cpp:122
unsigned nRows
The current number of rows, columns, and reserved columns.
Definition: Matrix.h:240
void removeColumn(unsigned pos)
Definition: Matrix.cpp:196
unsigned appendExtraRow()
Add an extra row at the bottom of the matrix and return its position.
Definition: Matrix.cpp:67
unsigned nReservedColumns
Definition: Matrix.h:240
void addToColumn(unsigned sourceColumn, unsigned targetColumn, const T &scale)
Add scale multiples of the source column to the target column.
Definition: Matrix.cpp:321
Matrix< T > getSubMatrix(unsigned fromRow, unsigned toRow, unsigned fromColumn, unsigned toColumn) const
Definition: Matrix.cpp:386
void print(raw_ostream &os) const
Print the matrix.
Definition: Matrix.cpp:402
void copyRow(unsigned sourceRow, unsigned targetRow)
Definition: Matrix.cpp:246
void scaleRow(unsigned row, const T &scale)
Multiply the specified row by a factor of scale.
Definition: Matrix.cpp:315
void insertColumn(unsigned pos)
Definition: Matrix.cpp:150
MutableArrayRef< T > getRow(unsigned row)
Get a [Mutable]ArrayRef corresponding to the specified row.
Definition: Matrix.cpp:132
void removeColumns(unsigned pos, unsigned count)
Remove the columns having positions pos, pos + 1, ...
Definition: Matrix.cpp:200
void insertColumns(unsigned pos, unsigned count)
Insert columns having positions pos, pos + 1, ...
Definition: Matrix.cpp:154
void setRow(unsigned row, ArrayRef< T > elems)
Set the specified row to elems.
Definition: Matrix.cpp:142
std::pair< Matrix< T >, Matrix< T > > splitByBitset(ArrayRef< int > indicator)
Split the rows of a matrix into two matrices according to which bits are 1 and which are 0 in a given...
Definition: Matrix.cpp:414
void removeRow(unsigned pos)
Definition: Matrix.cpp:232
bool operator==(const Matrix< T > &m) const
We cannot use the default implementation of operator== as it compares fields like reservedColumns etc...
Definition: Matrix.cpp:35
SmallVector< T, 16 > data
Stores the data.
Definition: Matrix.h:244
unsigned getNumColumns() const
Definition: Matrix.h:87
void resizeVertically(unsigned newNRows)
Definition: Matrix.cpp:106
unsigned getNumReservedRows() const
Return the maximum number of rows/columns that can be added without incurring a reallocation.
Definition: Matrix.cpp:57
Matrix< T > transpose() const
Definition: Matrix.cpp:82
SmallVector< T, 8 > preMultiplyWithRow(ArrayRef< T > rowVec) const
The given vector is interpreted as a row vector v.
Definition: Matrix.cpp:348
static Matrix identity(unsigned dimension)
Return the identity matrix of the specified dimension.
Definition: Matrix.cpp:49
void insertRow(unsigned pos)
Definition: Matrix.cpp:214
SmallVector< T, 8 > postMultiplyWithColumn(ArrayRef< T > colVec) const
The given vector is interpreted as a column vector v.
Definition: Matrix.cpp:359
void negateMatrix()
Negate the entire matrix.
Definition: Matrix.cpp:342
void swapRows(unsigned row, unsigned otherRow)
Swap the given rows.
Definition: Matrix.cpp:112
void resizeHorizontally(unsigned newNColumns)
Definition: Matrix.cpp:92
void reserveRows(unsigned rows)
Reserve enough space to resize to the specified number of rows without reallocations.
Definition: Matrix.cpp:62
void negateColumn(unsigned column)
Negate the specified column.
Definition: Matrix.cpp:330
void resize(unsigned newNRows, unsigned newNColumns)
Resize the matrix to the specified dimensions.
Definition: Matrix.cpp:100
void fillRow(unsigned row, const T &value)
Definition: Matrix.cpp:254
void addToRow(unsigned sourceRow, unsigned targetRow, const T &scale)
Add scale multiples of the source row to the target row.
Definition: Matrix.cpp:301
void negateRow(unsigned row)
Negate the specified row.
Definition: Matrix.cpp:336
T & at(unsigned row, unsigned column)
Access the element at the specified row and column.
Definition: Matrix.h:61
void removeRows(unsigned pos, unsigned count)
Remove the rows having positions pos, pos + 1, ...
Definition: Matrix.cpp:236
MPInt round(const Fraction &f)
Definition: Fraction.h:133
Fraction dotProduct(ArrayRef< Fraction > a, ArrayRef< Fraction > b)
Compute the dot product of two vectors.
Definition: Utils.cpp:533
MPInt normalizeRange(MutableArrayRef< MPInt > range)
Divide the range by its gcd and return the gcd.
Definition: Utils.cpp:356
LLVM_ATTRIBUTE_ALWAYS_INLINE MPInt lcm(const MPInt &a, const MPInt &b)
Returns the least common multiple of 'a' and 'b'.
Definition: MPInt.h:407
LLVM_ATTRIBUTE_ALWAYS_INLINE MPInt floorDiv(const MPInt &lhs, const MPInt &rhs)
Definition: MPInt.h:382
Include the generated interface declarations.
A class to represent fractions.
Definition: Fraction.h:28
MPInt getAsInteger() const
Definition: Fraction.h:48
Eliminates variable at the specified position using Fourier-Motzkin variable elimination.