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 | |
18 | using namespace mlir; |
19 | using namespace presburger; |
20 | |
21 | template <typename T> |
22 | Matrix<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. |
32 | template <typename T> |
33 | bool 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 | |
46 | template <typename T> |
47 | Matrix<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 | |
54 | template <typename T> |
55 | unsigned Matrix<T>::getNumReservedRows() const { |
56 | return data.capacity() / nReservedColumns; |
57 | } |
58 | |
59 | template <typename T> |
60 | void Matrix<T>::reserveRows(unsigned rows) { |
61 | data.reserve(rows * nReservedColumns); |
62 | } |
63 | |
64 | template <typename T> |
65 | unsigned Matrix<T>::appendExtraRow() { |
66 | resizeVertically(newNRows: nRows + 1); |
67 | return nRows - 1; |
68 | } |
69 | |
70 | template <typename T> |
71 | unsigned 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 | |
79 | template <typename T> |
80 | Matrix<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 | |
89 | template <typename T> |
90 | void 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 | |
97 | template <typename T> |
98 | void Matrix<T>::resize(unsigned newNRows, unsigned newNColumns) { |
99 | resizeHorizontally(newNColumns); |
100 | resizeVertically(newNRows); |
101 | } |
102 | |
103 | template <typename T> |
104 | void Matrix<T>::resizeVertically(unsigned newNRows) { |
105 | nRows = newNRows; |
106 | data.resize(nRows * nReservedColumns); |
107 | } |
108 | |
109 | template <typename T> |
110 | void 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 | |
119 | template <typename T> |
120 | void 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 | |
129 | template <typename T> |
130 | MutableArrayRef<T> Matrix<T>::getRow(unsigned row) { |
131 | return {&data[row * nReservedColumns], nColumns}; |
132 | } |
133 | |
134 | template <typename T> |
135 | ArrayRef<T> Matrix<T>::getRow(unsigned row) const { |
136 | return {&data[row * nReservedColumns], nColumns}; |
137 | } |
138 | |
139 | template <typename T> |
140 | void 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 | |
147 | template <typename T> |
148 | void Matrix<T>::insertColumn(unsigned pos) { |
149 | insertColumns(pos, count: 1); |
150 | } |
151 | template <typename T> |
152 | void 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 | |
193 | template <typename T> |
194 | void Matrix<T>::removeColumn(unsigned pos) { |
195 | removeColumns(pos, count: 1); |
196 | } |
197 | template <typename T> |
198 | void 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 | |
211 | template <typename T> |
212 | void Matrix<T>::insertRow(unsigned pos) { |
213 | insertRows(pos, count: 1); |
214 | } |
215 | template <typename T> |
216 | void 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 | |
229 | template <typename T> |
230 | void Matrix<T>::removeRow(unsigned pos) { |
231 | removeRows(pos, count: 1); |
232 | } |
233 | template <typename T> |
234 | void 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 | |
243 | template <typename T> |
244 | void 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 | |
251 | template <typename T> |
252 | void 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. |
265 | template <typename T> |
266 | void 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 | |
298 | template <typename T> |
299 | void Matrix<T>::addToRow(unsigned sourceRow, unsigned targetRow, |
300 | const T &scale) { |
301 | addToRow(targetRow, getRow(sourceRow), scale); |
302 | } |
303 | |
304 | template <typename T> |
305 | void 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 | |
312 | template <typename T> |
313 | void Matrix<T>::scaleRow(unsigned row, const T &scale) { |
314 | for (unsigned col = 0; col < nColumns; ++col) |
315 | at(row, col) *= scale; |
316 | } |
317 | |
318 | template <typename T> |
319 | void 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 | |
327 | template <typename T> |
328 | void Matrix<T>::negateColumn(unsigned column) { |
329 | for (unsigned row = 0, e = getNumRows(); row < e; ++row) |
330 | at(row, column) = -at(row, column); |
331 | } |
332 | |
333 | template <typename T> |
334 | void Matrix<T>::negateRow(unsigned row) { |
335 | for (unsigned column = 0, e = getNumColumns(); column < e; ++column) |
336 | at(row, column) = -at(row, column); |
337 | } |
338 | |
339 | template <typename T> |
340 | void Matrix<T>::negateMatrix() { |
341 | for (unsigned row = 0; row < nRows; ++row) |
342 | negateRow(row); |
343 | } |
344 | |
345 | template <typename T> |
346 | SmallVector<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 | |
356 | template <typename T> |
357 | SmallVector<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. |
373 | static 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 | |
383 | template <typename T> |
384 | Matrix<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 | |
399 | template <typename T> |
400 | void 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. |
416 | template <typename T> |
417 | std::pair<Matrix<T>, Matrix<T>> |
418 | Matrix<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 | |
429 | template <typename T> |
430 | void Matrix<T>::dump() const { |
431 | print(os&: llvm::errs()); |
432 | } |
433 | |
434 | template <typename T> |
435 | bool 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 | |
449 | namespace mlir { |
450 | namespace presburger { |
451 | template class Matrix<DynamicAPInt>; |
452 | template class Matrix<Fraction>; |
453 | } // namespace presburger |
454 | } // namespace mlir |
455 | |
456 | IntMatrix 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 | |
463 | std::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 | |
549 | DynamicAPInt IntMatrix::normalizeRow(unsigned row, unsigned cols) { |
550 | return normalizeRange(range: getRow(row).slice(N: 0, M: cols)); |
551 | } |
552 | |
553 | DynamicAPInt IntMatrix::normalizeRow(unsigned row) { |
554 | return normalizeRow(row, cols: getNumColumns()); |
555 | } |
556 | |
557 | DynamicAPInt 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 | |
580 | FracMatrix FracMatrix::identity(unsigned dimension) { |
581 | return Matrix::identity(dimension); |
582 | } |
583 | |
584 | FracMatrix::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 | |
591 | Fraction 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 | |
676 | FracMatrix 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. |
724 | void 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 | |
764 | IntMatrix 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 |
Definitions
- Matrix
- operator==
- identity
- getNumReservedRows
- reserveRows
- appendExtraRow
- appendExtraRow
- transpose
- resizeHorizontally
- resize
- resizeVertically
- swapRows
- swapColumns
- getRow
- getRow
- setRow
- insertColumn
- insertColumns
- removeColumn
- removeColumns
- insertRow
- insertRows
- removeRow
- removeRows
- copyRow
- fillRow
- moveColumns
- addToRow
- addToRow
- scaleRow
- addToColumn
- negateColumn
- negateRow
- negateMatrix
- preMultiplyWithRow
- postMultiplyWithColumn
- modEntryColumnOperation
- getSubMatrix
- splitByBitset
- dump
- hasConsistentState
- identity
- computeHermiteNormalForm
- normalizeRow
- normalizeRow
- determinant
- identity
- FracMatrix
- determinant
- gramSchmidt
- LLL
Improve your Profiling and Debugging skills
Find out more