12 #include "llvm/Support/MathExtras.h"
13 #include "llvm/Support/raw_ostream.h"
19 using namespace presburger;
23 unsigned reservedColumns)
24 : nRows(
rows), nColumns(columns),
25 nReservedColumns(std::
max(nColumns, reservedColumns)),
26 data(nRows * nReservedColumns) {
34 if (nRows != m.getNumRows())
36 if (nColumns != m.getNumColumns())
39 for (
unsigned i = 0; i < nRows; i++)
40 if (getRow(i) != m.getRow(i))
48 Matrix matrix(dimension, dimension);
49 for (
unsigned i = 0; i < dimension; ++i)
56 return data.capacity() / nReservedColumns;
61 data.reserve(
rows * nReservedColumns);
66 resizeVertically(nRows + 1);
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];
82 for (
unsigned row = 0; row < nRows; ++row)
83 for (
unsigned col = 0; col < nColumns; ++col)
84 transp(col, row) = at(row, col);
91 if (newNColumns < nColumns)
92 removeColumns(newNColumns, nColumns - newNColumns);
93 if (newNColumns > nColumns)
94 insertColumns(nColumns, newNColumns - nColumns);
99 resizeHorizontally(newNColumns);
100 resizeVertically(newNRows);
103 template <
typename T>
106 data.resize(nRows * nReservedColumns);
109 template <
typename T>
111 assert((row < getNumRows() && otherRow < getNumRows()) &&
112 "Given row out of bounds");
115 for (
unsigned col = 0; col < nColumns; col++)
116 std::swap(at(row, col), at(otherRow, col));
119 template <
typename T>
121 assert((column < getNumColumns() && otherColumn < getNumColumns()) &&
122 "Given column out of bounds");
123 if (column == otherColumn)
125 for (
unsigned row = 0; row < nRows; row++)
126 std::swap(at(row, column), at(row, otherColumn));
129 template <
typename T>
131 return {&data[row * nReservedColumns], nColumns};
134 template <
typename T>
136 return {&data[row * nReservedColumns], nColumns};
139 template <
typename T>
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];
147 template <
typename T>
149 insertColumns(pos, 1);
151 template <
typename T>
155 assert(pos <= nColumns);
156 unsigned oldNReservedColumns = nReservedColumns;
157 if (nColumns + count > nReservedColumns) {
158 nReservedColumns = llvm::NextPowerOf2(nColumns + count);
159 data.resize(nRows * nReservedColumns);
163 for (
int ri = nRows - 1; ri >= 0; --ri) {
164 for (
int ci = nReservedColumns - 1; ci >= 0; --ci) {
167 T &dest = data[r * nReservedColumns + c];
175 }
else if (c >= pos + count) {
177 dest = data[r * oldNReservedColumns + c - count];
178 }
else if (c >= pos) {
185 if (nReservedColumns == oldNReservedColumns)
187 dest = data[r * oldNReservedColumns + c];
193 template <
typename T>
195 removeColumns(pos, 1);
197 template <
typename T>
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)
211 template <
typename T>
215 template <
typename T>
220 assert(pos <= nRows);
221 resizeVertically(nRows + count);
222 for (
int r = nRows - 1; r >= int(pos + count); --r)
223 copyRow(r - count, r);
224 for (
int r = pos + count - 1; r >= int(pos); --r)
225 for (
unsigned c = 0; c < nColumns; ++c)
229 template <
typename T>
233 template <
typename T>
237 assert(pos + count - 1 <= nRows);
238 for (
unsigned r = pos; r + count < nRows; ++r)
239 copyRow(r + count, r);
240 resizeVertically(nRows - count);
243 template <
typename T>
245 if (sourceRow == targetRow)
247 for (
unsigned c = 0; c < nColumns; ++c)
248 at(targetRow, c) = at(sourceRow, c);
251 template <
typename T>
253 for (
unsigned col = 0; col < nColumns; ++col)
254 at(row, col) = value;
265 template <
typename T>
270 int offset = dstPos - srcPos;
274 assert(srcPos + num <= getNumColumns() &&
275 "move source range exceeds matrix columns");
276 assert(dstPos + num <= getNumColumns() &&
277 "move destination range exceeds matrix columns");
279 unsigned insertCount = offset > 0 ? offset : -offset;
280 unsigned finalAdjStart = offset > 0 ? srcPos : srcPos + num;
281 unsigned curAdjStart = offset > 0 ? srcPos + num : dstPos;
285 insertColumns(finalAdjStart, insertCount);
287 if (finalAdjStart < curAdjStart)
288 curAdjStart += insertCount;
291 for (
unsigned i = 0; i < insertCount; ++i)
292 swapColumns(finalAdjStart + i, curAdjStart + i);
295 removeColumns(curAdjStart, insertCount);
298 template <
typename T>
301 addToRow(targetRow, getRow(sourceRow), scale);
304 template <
typename T>
308 for (
unsigned col = 0; col < nColumns; ++col)
309 at(row, col) += scale * rowVec[col];
312 template <
typename T>
314 for (
unsigned col = 0; col < nColumns; ++col)
315 at(row, col) *= scale;
318 template <
typename T>
323 for (
unsigned row = 0, e = getNumRows(); row < e; ++row)
324 at(row, targetColumn) += scale * at(row, sourceColumn);
327 template <
typename T>
329 for (
unsigned row = 0, e = getNumRows(); row < e; ++row)
330 at(row, column) = -at(row, column);
333 template <
typename T>
335 for (
unsigned column = 0, e = getNumColumns(); column < e; ++column)
336 at(row, column) = -at(row, column);
339 template <
typename T>
341 for (
unsigned row = 0; row < nRows; ++row)
345 template <
typename T>
347 assert(rowVec.size() == getNumRows() &&
"Invalid row vector dimension!");
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);
356 template <
typename T>
358 assert(getNumColumns() == colVec.size() &&
359 "Invalid column vector dimension!");
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];
374 unsigned sourceCol,
unsigned targetCol,
376 assert(m(row, sourceCol) != 0 &&
"Cannot divide by zero!");
377 assert(m(row, sourceCol) > 0 &&
"Source must be positive!");
378 DynamicAPInt ratio = -floorDiv(m(row, targetCol), m(row, sourceCol));
379 m.addToColumn(sourceCol, targetCol, ratio);
380 otherMatrix.
addToColumn(sourceCol, targetCol, ratio);
383 template <
typename T>
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);
399 template <
typename T>
401 for (
unsigned row = 0; row < nRows; ++row) {
402 for (
unsigned column = 0; column < nColumns; ++column)
403 os << at(row, column) <<
' ';
410 template <
typename T>
413 Matrix<T> rowsForOne(0, nColumns), rowsForZero(0, nColumns);
414 for (
unsigned i = 0; i < nRows; i++) {
415 if (indicator[i] == 1)
416 rowsForOne.appendExtraRow(getRow(i));
420 return {rowsForOne, rowsForZero};
423 template <
typename T>
428 template <
typename T>
430 if (data.size() != nRows * nReservedColumns)
432 if (nColumns > nReservedColumns)
434 #ifdef EXPENSIVE_CHECKS
435 for (
unsigned r = 0; r < nRows; ++r)
436 for (
unsigned c = nColumns; c < nReservedColumns; ++c)
437 if (data[r * nReservedColumns + c] != 0)
444 namespace presburger {
452 for (
unsigned i = 0; i < dimension; ++i)
464 unsigned echelonCol = 0;
469 for (
unsigned row = 0; row < h.
getNumRows(); ++row) {
471 unsigned nonZeroCol = echelonCol;
472 for (
unsigned e = h.
getNumColumns(); nonZeroCol < e; ++nonZeroCol) {
473 if (h(row, nonZeroCol) == 0)
485 if (nonZeroCol != echelonCol) {
491 if (h(row, echelonCol) < 0) {
497 for (
unsigned i = echelonCol + 1, e = h.
getNumColumns(); i < e; ++i) {
507 unsigned targetCol = i, sourceCol = echelonCol;
519 while (h(row, targetCol) != 0 && h(row, sourceCol) != 0) {
521 std::swap(targetCol, sourceCol);
526 if (h(row, echelonCol) == 0) {
534 for (
unsigned i = 0; i < echelonCol; ++i)
553 "determinant can only be calculated for square matrices!");
558 DynamicAPInt detM = m.determinant(&fracInverse).getAsInteger();
561 return DynamicAPInt(0);
567 for (
unsigned i = 0; i <
nRows; i++)
569 inverse->
at(i,
j) = (fracInverse.
at(i,
j) * detM).getAsInteger();
579 :
FracMatrix(m.getNumRows(), m.getNumColumns()) {
580 for (
unsigned i = 0, r = m.getNumRows(); i < r; i++)
581 for (
unsigned j = 0, c = m.getNumColumns();
j < c;
j++)
582 this->
at(i,
j) = m.at(i,
j);
587 "determinant can only be calculated for square matrices!");
603 for (
unsigned i = 0; i <
nRows; i++) {
608 for (
unsigned j = i + 1;
j <
nRows;
j++) {
624 for (
unsigned j = 0;
j < i;
j++) {
631 m.addToRow(i,
j, -a / b);
638 for (
unsigned j = i + 1;
j <
nRows;
j++) {
645 m.addToRow(i,
j, -a / b);
656 for (
unsigned i = 0; i <
nRows; i++)
657 for (
unsigned j = 0;
j <
nRows;
j++)
658 tempInv.
at(i,
j) = tempInv.
at(i,
j) / m(i, i);
660 *inverse = std::move(tempInv);
664 for (
unsigned i = 0; i <
nRows; i++)
679 for (
unsigned i = 1, e =
getNumRows(); i < e; i++) {
680 for (
unsigned j = 0;
j < i;
j++) {
682 assert(jNormSquared != 0 &&
"some row became zero! Inputs to this "
683 "function must be linearly independent.");
719 DynamicAPInt nearest;
731 for (
unsigned j = k - 1;
j < k;
j--) {
753 k = k > 1 ? k - 1 : 1;
761 IntMatrix normalized(numRows, numColumns);
763 DynamicAPInt lcmDenoms = DynamicAPInt(1);
764 for (
unsigned i = 0; i < numRows; i++) {
766 for (
unsigned j = 0;
j < numColumns;
j++)
767 lcmDenoms = lcm(lcmDenoms,
at(i,
j).den);
769 for (
unsigned j = 0;
j < numColumns;
j++)
770 normalized(i,
j) = (
at(i,
j) * lcmDenoms).getAsInteger();
static void modEntryColumnOperation(Matrix< DynamicAPInt > &m, unsigned row, unsigned sourceCol, unsigned targetCol, Matrix< DynamicAPInt > &otherMatrix)
Set M(row, targetCol) to its remainder on division by M(row, sourceCol) by subtracting from column ta...
static Value max(ImplicitLocOpBuilder &builder, Value value, Value bound)
static void print(spirv::VerCapExtAttr triple, DialectAsmPrinter &printer)
FracMatrix gramSchmidt() const
IntMatrix normalizeRows() const
static FracMatrix identity(unsigned dimension)
Return the identity matrix of the specified dimension.
FracMatrix(unsigned rows, unsigned columns, unsigned reservedRows=0, unsigned reservedColumns=0)
Fraction determinant(FracMatrix *inverse=nullptr) const
static IntMatrix identity(unsigned dimension)
Return the identity matrix of the specified dimension.
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...
DynamicAPInt determinant(IntMatrix *inverse=nullptr) const
IntMatrix(unsigned rows, unsigned columns, unsigned reservedRows=0, unsigned reservedColumns=0)
DynamicAPInt normalizeRow(unsigned row, unsigned nCols)
Divide the first nCols of the specified row by their GCD.
This is a class to represent a resizable matrix.
void moveColumns(unsigned srcPos, unsigned num, unsigned dstPos)
Move the columns in the source range [srcPos, srcPos + num) to the specified destination [dstPos,...
bool hasConsistentState() const
Return whether the Matrix is in a consistent state with all its invariants satisfied.
void insertRows(unsigned pos, unsigned count)
Insert rows having positions pos, pos + 1, ...
unsigned getNumRows() const
void swapColumns(unsigned column, unsigned otherColumn)
Swap the given columns.
unsigned nRows
The current number of rows, columns, and reserved columns.
void removeColumn(unsigned pos)
unsigned appendExtraRow()
Add an extra row at the bottom of the matrix and return its position.
unsigned nReservedColumns
void addToColumn(unsigned sourceColumn, unsigned targetColumn, const T &scale)
Add scale multiples of the source column to the target column.
Matrix< T > getSubMatrix(unsigned fromRow, unsigned toRow, unsigned fromColumn, unsigned toColumn) const
void print(raw_ostream &os) const
Print the matrix.
void copyRow(unsigned sourceRow, unsigned targetRow)
void scaleRow(unsigned row, const T &scale)
Multiply the specified row by a factor of scale.
void insertColumn(unsigned pos)
MutableArrayRef< T > getRow(unsigned row)
Get a [Mutable]ArrayRef corresponding to the specified row.
void removeColumns(unsigned pos, unsigned count)
Remove the columns having positions pos, pos + 1, ...
void insertColumns(unsigned pos, unsigned count)
Insert columns having positions pos, pos + 1, ...
void setRow(unsigned row, ArrayRef< T > elems)
Set the specified row to elems.
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...
void removeRow(unsigned pos)
bool operator==(const Matrix< T > &m) const
We cannot use the default implementation of operator== as it compares fields like reservedColumns etc...
SmallVector< T, 16 > data
Stores the data.
unsigned getNumColumns() const
void resizeVertically(unsigned newNRows)
unsigned getNumReservedRows() const
Return the maximum number of rows/columns that can be added without incurring a reallocation.
Matrix< T > transpose() const
SmallVector< T, 8 > preMultiplyWithRow(ArrayRef< T > rowVec) const
The given vector is interpreted as a row vector v.
static Matrix identity(unsigned dimension)
Return the identity matrix of the specified dimension.
void insertRow(unsigned pos)
SmallVector< T, 8 > postMultiplyWithColumn(ArrayRef< T > colVec) const
The given vector is interpreted as a column vector v.
void negateMatrix()
Negate the entire matrix.
void swapRows(unsigned row, unsigned otherRow)
Swap the given rows.
void resizeHorizontally(unsigned newNColumns)
void reserveRows(unsigned rows)
Reserve enough space to resize to the specified number of rows without reallocations.
void negateColumn(unsigned column)
Negate the specified column.
void resize(unsigned newNRows, unsigned newNColumns)
Resize the matrix to the specified dimensions.
void fillRow(unsigned row, const T &value)
void addToRow(unsigned sourceRow, unsigned targetRow, const T &scale)
Add scale multiples of the source row to the target row.
void negateRow(unsigned row)
Negate the specified row.
T & at(unsigned row, unsigned column)
Access the element at the specified row and column.
void removeRows(unsigned pos, unsigned count)
Remove the rows having positions pos, pos + 1, ...
DynamicAPInt round(const Fraction &f)
Fraction dotProduct(ArrayRef< Fraction > a, ArrayRef< Fraction > b)
Compute the dot product of two vectors.
DynamicAPInt normalizeRange(MutableArrayRef< DynamicAPInt > range)
Divide the range by its gcd and return the gcd.
Include the generated interface declarations.
A class to represent fractions.
Eliminates variable at the specified position using Fourier-Motzkin variable elimination.