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

source code of mlir/lib/Analysis/Presburger/Matrix.cpp