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/MPInt.h"
12#include "mlir/Analysis/Presburger/Utils.h"
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
20using namespace mlir;
21using namespace presburger;
22
23template <typename T>
24Matrix<T>::Matrix(unsigned rows, unsigned columns, unsigned reservedRows,
25 unsigned reservedColumns)
26 : nRows(rows), nColumns(columns),
27 nReservedColumns(std::max(a: nColumns, b: reservedColumns)),
28 data(nRows * nReservedColumns) {
29 data.reserve(std::max(a: nRows, b: 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.
34template <typename T>
35bool 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
48template <typename T>
49Matrix<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
56template <typename T>
57unsigned Matrix<T>::getNumReservedRows() const {
58 return data.capacity() / nReservedColumns;
59}
60
61template <typename T>
62void Matrix<T>::reserveRows(unsigned rows) {
63 data.reserve(rows * nReservedColumns);
64}
65
66template <typename T>
67unsigned Matrix<T>::appendExtraRow() {
68 resizeVertically(newNRows: nRows + 1);
69 return nRows - 1;
70}
71
72template <typename T>
73unsigned Matrix<T>::appendExtraRow(ArrayRef<T> elems) {
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
81template <typename T>
82Matrix<T> Matrix<T>::transpose() const {
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
91template <typename T>
92void Matrix<T>::resizeHorizontally(unsigned newNColumns) {
93 if (newNColumns < nColumns)
94 removeColumns(pos: newNColumns, count: nColumns - newNColumns);
95 if (newNColumns > nColumns)
96 insertColumns(pos: nColumns, count: newNColumns - nColumns);
97}
98
99template <typename T>
100void Matrix<T>::resize(unsigned newNRows, unsigned newNColumns) {
101 resizeHorizontally(newNColumns);
102 resizeVertically(newNRows);
103}
104
105template <typename T>
106void Matrix<T>::resizeVertically(unsigned newNRows) {
107 nRows = newNRows;
108 data.resize(nRows * nReservedColumns);
109}
110
111template <typename T>
112void 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
121template <typename T>
122void 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
131template <typename T>
132MutableArrayRef<T> Matrix<T>::getRow(unsigned row) {
133 return {&data[row * nReservedColumns], nColumns};
134}
135
136template <typename T>
137ArrayRef<T> Matrix<T>::getRow(unsigned row) const {
138 return {&data[row * nReservedColumns], nColumns};
139}
140
141template <typename T>
142void 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
149template <typename T>
150void Matrix<T>::insertColumn(unsigned pos) {
151 insertColumns(pos, count: 1);
152}
153template <typename T>
154void 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(A: 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
195template <typename T>
196void Matrix<T>::removeColumn(unsigned pos) {
197 removeColumns(pos, count: 1);
198}
199template <typename T>
200void 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
213template <typename T>
214void Matrix<T>::insertRow(unsigned pos) {
215 insertRows(pos, count: 1);
216}
217template <typename T>
218void Matrix<T>::insertRows(unsigned pos, unsigned count) {
219 if (count == 0)
220 return;
221
222 assert(pos <= nRows);
223 resizeVertically(newNRows: nRows + count);
224 for (int r = nRows - 1; r >= int(pos + count); --r)
225 copyRow(sourceRow: r - count, targetRow: 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
231template <typename T>
232void Matrix<T>::removeRow(unsigned pos) {
233 removeRows(pos, count: 1);
234}
235template <typename T>
236void 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(sourceRow: r + count, targetRow: r);
242 resizeVertically(newNRows: nRows - count);
243}
244
245template <typename T>
246void 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
253template <typename T>
254void 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.
267template <typename T>
268void 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(pos: finalAdjStart, count: 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(column: finalAdjStart + i, otherColumn: curAdjStart + i);
295
296 // Delete the now redundant zero columns.
297 removeColumns(pos: curAdjStart, count: insertCount);
298}
299
300template <typename T>
301void Matrix<T>::addToRow(unsigned sourceRow, unsigned targetRow,
302 const T &scale) {
303 addToRow(targetRow, getRow(sourceRow), scale);
304}
305
306template <typename T>
307void 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
314template <typename T>
315void Matrix<T>::scaleRow(unsigned row, const T &scale) {
316 for (unsigned col = 0; col < nColumns; ++col)
317 at(row, col) *= scale;
318}
319
320template <typename T>
321void 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
329template <typename T>
330void Matrix<T>::negateColumn(unsigned column) {
331 for (unsigned row = 0, e = getNumRows(); row < e; ++row)
332 at(row, column) = -at(row, column);
333}
334
335template <typename T>
336void Matrix<T>::negateRow(unsigned row) {
337 for (unsigned column = 0, e = getNumColumns(); column < e; ++column)
338 at(row, column) = -at(row, column);
339}
340
341template <typename T>
342void Matrix<T>::negateMatrix() {
343 for (unsigned row = 0; row < nRows; ++row)
344 negateRow(row);
345}
346
347template <typename T>
348SmallVector<T, 8> Matrix<T>::preMultiplyWithRow(ArrayRef<T> rowVec) const {
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
358template <typename T>
359SmallVector<T, 8> Matrix<T>::postMultiplyWithColumn(ArrayRef<T> colVec) const {
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.
375static 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(lhs: m(row, targetCol), rhs: m(row, sourceCol));
381 m.addToColumn(sourceColumn: sourceCol, targetColumn: targetCol, scale: ratio);
382 otherMatrix.addToColumn(sourceColumn: sourceCol, targetColumn: targetCol, scale: ratio);
383}
384
385template <typename T>
386Matrix<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
401template <typename T>
402void 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.
412template <typename T>
413std::pair<Matrix<T>, Matrix<T>>
414Matrix<T>::splitByBitset(ArrayRef<int> indicator) {
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
425template <typename T>
426void Matrix<T>::dump() const {
427 print(os&: llvm::errs());
428}
429
430template <typename T>
431bool Matrix<T>::hasConsistentState() const {
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
445namespace mlir {
446namespace presburger {
447template class Matrix<MPInt>;
448template class Matrix<Fraction>;
449} // namespace presburger
450} // namespace mlir
451
452IntMatrix 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
459std::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;
464 IntMatrix u = IntMatrix::identity(dimension: h.getNumColumns());
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(column: nonZeroCol, otherColumn: echelonCol);
489 u.swapColumns(column: nonZeroCol, otherColumn: echelonCol);
490 }
491
492 // Make h(row, echelonCol) non-negative.
493 if (h(row, echelonCol) < 0) {
494 h.negateColumn(column: echelonCol);
495 u.negateColumn(column: 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(column: i);
506 u.negateColumn(column: 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(m&: h, row, sourceCol, targetCol, otherMatrix&: u);
523 std::swap(a&: targetCol, b&: 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(column: i, otherColumn: echelonCol);
530 u.swapColumns(column: i, otherColumn: 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(m&: h, row, sourceCol: echelonCol, targetCol: i, otherMatrix&: u);
538
539 ++echelonCol;
540 }
541
542 return {h, u};
543}
544
545MPInt IntMatrix::normalizeRow(unsigned row, unsigned cols) {
546 return normalizeRange(range: getRow(row).slice(N: 0, M: cols));
547}
548
549MPInt IntMatrix::normalizeRow(unsigned row) {
550 return normalizeRow(row, cols: getNumColumns());
551}
552
553MPInt IntMatrix::determinant(IntMatrix *inverse) const {
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(inverse: &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(row: i, column: j) = (fracInverse.at(row: i, column: j) * detM).getAsInteger();
572
573 return detM;
574}
575
576FracMatrix FracMatrix::identity(unsigned dimension) {
577 return Matrix::identity(dimension);
578}
579
580FracMatrix::FracMatrix(IntMatrix m)
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(row: i, column: j) = m.at(row: i, column: j);
585}
586
587Fraction FracMatrix::determinant(FracMatrix *inverse) const {
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(dimension: 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(row: j, otherRow: i);
613 if (inverse)
614 tempInv.swapRows(row: j, otherRow: i);
615 break;
616 }
617 }
618
619 b = m.at(row: i, column: 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(row: j, column: i) == 0)
628 continue;
629 a = m.at(row: j, column: i);
630 // Set element (j, i) to zero
631 // by subtracting the ith row,
632 // appropriately scaled.
633 m.addToRow(sourceRow: i, targetRow: j, scale: -a / b);
634 tempInv.addToRow(sourceRow: i, targetRow: j, scale: -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(row: j, column: i) == 0)
642 continue;
643 a = m.at(row: j, column: i);
644 // Set element (j, i) to zero
645 // by subtracting the ith row,
646 // appropriately scaled.
647 m.addToRow(sourceRow: i, targetRow: j, scale: -a / b);
648 if (inverse)
649 tempInv.addToRow(sourceRow: i, targetRow: j, scale: -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(row: i, column: j) = tempInv.at(row: i, column: 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(row: i, column: i);
668
669 return determinant;
670}
671
672FracMatrix FracMatrix::gramSchmidt() const {
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(a: orth.getRow(row: j), b: orth.getRow(row: j));
684 assert(jNormSquared != 0 && "some row became zero! Inputs to this "
685 "function must be linearly independent.");
686 Fraction projectionScale =
687 dotProduct(a: orth.getRow(row: i), b: orth.getRow(row: j)) / jNormSquared;
688 orth.addToRow(sourceRow: j, targetRow: i, scale: -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.
720void FracMatrix::LLL(Fraction delta) {
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(a: getRow(row: k), b: gsOrth.getRow(row: j)) /
736 dotProduct(a: gsOrth.getRow(row: j), b: gsOrth.getRow(row: j));
737 nearest = round(f: mu);
738 // Subtract b_j scaled by the integer nearest to μ_jk from b_k.
739 addToRow(row: k, rowVec: getRow(row: j), scale: -Fraction(nearest, 1));
740 gsOrth = gramSchmidt(); // Update orthogonalization.
741 }
742 mu = dotProduct(a: getRow(row: k), b: gsOrth.getRow(row: k - 1)) /
743 dotProduct(a: gsOrth.getRow(row: k - 1), b: gsOrth.getRow(row: k - 1));
744 // Check the Lovasz condition for b_k and b_{k-1}.
745 if (dotProduct(a: gsOrth.getRow(row: k), b: gsOrth.getRow(row: k)) >
746 (delta - mu * mu) *
747 dotProduct(a: gsOrth.getRow(row: k - 1), b: gsOrth.getRow(row: 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(row: k, otherRow: k - 1);
754 gsOrth = gramSchmidt(); // Update orthogonalization.
755 k = k > 1 ? k - 1 : 1;
756 }
757 }
758}
759
760IntMatrix FracMatrix::normalizeRows() const {
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(a: lcmDenoms, b: at(row: i, column: 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(row: i, column: j) * lcmDenoms).getAsInteger();
773 }
774 return normalized;
775}
776

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