Matrix.cpp Source File (original) (raw)

12 #include "llvm/Support/MathExtras.h"

13 #include "llvm/Support/raw_ostream.h"

14 #include

15 #include

16 #include

18 using namespace mlir;

19 using namespace presburger;

21 template

23 unsigned reservedColumns)

24 : nRows(rows), nColumns(columns),

25 nReservedColumns(std::max(nColumns, reservedColumns)),

26 data(nRows * nReservedColumns) {

28 }

32 template

34 if (nRows != m.getNumRows())

35 return false;

36 if (nColumns != m.getNumColumns())

37 return false;

39 for (unsigned i = 0; i < nRows; i++)

40 if (getRow(i) != m.getRow(i))

41 return false;

43 return true;

44 }

46 template

48 Matrix matrix(dimension, dimension);

49 for (unsigned i = 0; i < dimension; ++i)

50 matrix(i, i) = 1;

51 return matrix;

52 }

54 template

56 return data.capacity() / nReservedColumns;

57 }

59 template

61 data.reserve(rows * nReservedColumns);

62 }

64 template

66 resizeVertically(nRows + 1);

67 return nRows - 1;

68 }

70 template

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 }

79 template

82 for (unsigned row = 0; row < nRows; ++row)

83 for (unsigned col = 0; col < nColumns; ++col)

84 transp(col, row) = at(row, col);

86 return transp;

87 }

89 template

91 if (newNColumns < nColumns)

92 removeColumns(newNColumns, nColumns - newNColumns);

93 if (newNColumns > nColumns)

94 insertColumns(nColumns, newNColumns - nColumns);

95 }

97 template

99 resizeHorizontally(newNColumns);

100 resizeVertically(newNRows);

102

103 template

105 nRows = newNRows;

106 data.resize(nRows * nReservedColumns);

107 }

108

109 template

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));

119 template

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

131 return {&data[row * nReservedColumns], nColumns};

132 }

134 template

136 return {&data[row * nReservedColumns], nColumns};

137 }

139 template

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

149 insertColumns(pos, 1);

151 template

153 if (count == 0)

154 return;

155 assert(pos <= nColumns);

156 unsigned oldNReservedColumns = nReservedColumns;

157 if (nColumns + count > nReservedColumns) {

158 nReservedColumns = llvm::NextPowerOf2(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) {

169

170

172

173

174 dest = 0;

175 } else if (c >= pos + count) {

176

177 dest = data[r * oldNReservedColumns + c - count];

178 } else if (c >= pos) {

179

180 dest = 0;

181 } else {

182

183

184

185 if (nReservedColumns == oldNReservedColumns)

187 dest = data[r * oldNReservedColumns + c];

188 }

190 }

191 }

193 template

195 removeColumns(pos, 1);

196 }

197 template

199 if (count == 0)

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

213 insertRows(pos, 1);

214 }

215 template

217 if (count == 0)

218 return;

219

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)

226 at(r, c) = 0;

227 }

228

229 template

231 removeRows(pos, 1);

232 }

233 template

235 if (count == 0)

236 return;

237 assert(pos + count - 1 <= nRows);

238 for (unsigned r = pos; r + count < nRows; ++r)

239 copyRow(r + count, r);

240 resizeVertically(nRows - count);

241 }

242

243 template

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

253 for (unsigned col = 0; col < nColumns; ++col)

254 at(row, col) = value;

255 }

256

257

258

259

260

261

262

263

264

265 template

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

283

284

285 insertColumns(finalAdjStart, insertCount);

286

287 if (finalAdjStart < curAdjStart)

288 curAdjStart += insertCount;

289

290

291 for (unsigned i = 0; i < insertCount; ++i)

292 swapColumns(finalAdjStart + i, curAdjStart + i);

293

294

295 removeColumns(curAdjStart, insertCount);

296 }

297

298 template

300 const T &scale) {

301 addToRow(targetRow, getRow(sourceRow), scale);

302 }

303

304 template

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

314 for (unsigned col = 0; col < nColumns; ++col)

315 at(row, col) *= scale;

316 }

317

318 template

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

329 for (unsigned row = 0, e = getNumRows(); row < e; ++row)

330 at(row, column) = -at(row, column);

331 }

332

333 template

335 for (unsigned column = 0, e = getNumColumns(); column < e; ++column)

336 at(row, column) = -at(row, column);

337 }

338

339 template

341 for (unsigned row = 0; row < nRows; ++row)

342 negateRow(row);

343 }

344

345 template

347 assert(rowVec.size() == getNumRows() && "Invalid row vector dimension!");

348

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

358 assert(getNumColumns() == colVec.size() &&

359 "Invalid column vector dimension!");

360

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

369

370

371

372

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);

381 }

382

383 template

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 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

402 for (unsigned row = 0; row < nRows; ++row)

403 for (unsigned column = 0; column < nColumns; ++column)

404 updatePrintMetrics(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(os, at(row, column), MIN_SPACING, ptm);

409 }

410 os << "\n";

411 }

412 }

413

414

415

416 template

419 Matrix 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

425 }

426 return {rowsForOne, rowsForZero};

427 }

428

429 template

431 print(llvm::errs());

432 }

433

434 template

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 {

453 }

454 }

455

457 IntMatrix matrix(dimension, dimension);

458 for (unsigned i = 0; i < dimension; ++i)

459 matrix(i, i) = 1;

460 return matrix;

461 }

462

464

465

466

469

470 unsigned echelonCol = 0;

471

472

473

474

475 for (unsigned row = 0; row < h.getNumRows(); ++row) {

476

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

485

487 continue;

488

489

490

491 if (nonZeroCol != echelonCol) {

494 }

495

496

497 if (h(row, echelonCol) < 0) {

500 }

501

502

503 for (unsigned i = echelonCol + 1, e = h.getNumColumns(); i < e; ++i) {

504

505

506

507

508 if (h(row, i) < 0) {

511 }

512

513 unsigned targetCol = i, sourceCol = echelonCol;

514

515

516

517

518

519

520

521

522

523

524

525 while (h(row, targetCol) != 0 && h(row, sourceCol) != 0) {

527 std::swap(targetCol, sourceCol);

528 }

529

530

531

532 if (h(row, echelonCol) == 0) {

535 }

536 }

537

538

539

540 for (unsigned i = 0; i < echelonCol; ++i)

542

543 ++echelonCol;

544 }

545

546 return {h, u};

547 }

548

551 }

552

555 }

556

559 "determinant can only be calculated for square matrices!");

560

562

564 DynamicAPInt detM = m.determinant(&fracInverse).getAsInteger();

565

566 if (detM == 0)

567 return DynamicAPInt(0);

568

569 if (!inverse)

570 return detM;

571

573 for (unsigned i = 0; i < nRows; i++)

575 inverse->at(i, j) = (fracInverse.at(i, j) * detM).getAsInteger();

576

577 return detM;

578 }

579

582 }

583

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(i, j) = m.at(i, j);

589 }

590

593 "determinant can only be calculated for square matrices!");

594

597 if (inverse)

599

601

602

603

604

605

606

607

608

609 for (unsigned i = 0; i < nRows; i++) {

610 if (m(i, i) == 0)

611

612

613

614 for (unsigned j = i + 1; j < nRows; j++) {

615 if (m(j, i) != 0) {

616 m.swapRows(j, i);

617 if (inverse)

619 break;

620 }

621 }

622

623 b = m.at(i, i);

624 if (b == 0)

625 return 0;

626

627

628

629 if (inverse) {

630 for (unsigned j = 0; j < i; j++) {

631 if (m.at(j, i) == 0)

632 continue;

633 a = m.at(j, i);

634

635

636

637 m.addToRow(i, j, -a / b);

639 }

640 }

641

642

643

644 for (unsigned j = i + 1; j < nRows; j++) {

645 if (m.at(j, i) == 0)

646 continue;

647 a = m.at(j, i);

648

649

650

651 m.addToRow(i, j, -a / b);

652 if (inverse)

654 }

655 }

656

657

658

659

660

661 if (inverse) {

662 for (unsigned i = 0; i < nRows; i++)

663 for (unsigned j = 0; j < nRows; j++)

664 tempInv.at(i, j) = tempInv.at(i, j) / m(i, i);

665

666 *inverse = std::move(tempInv);

667 }

668

670 for (unsigned i = 0; i < nRows; i++)

672

674 }

675

677

678

680

681

682

683

684

685 for (unsigned i = 1, e = getNumRows(); i < e; i++) {

686 for (unsigned j = 0; j < i; j++) {

688 assert(jNormSquared != 0 && "some row became zero! Inputs to this "

689 "function must be linearly independent.");

692 orth.addToRow(j, i, -projectionScale);

693 }

694 }

695 return orth;

696 }

697

698

699

700

701

702

703

704

705

706

707

708

709

710

711

712

713

714

715

716

717

718

719

720

721

722

723

725 DynamicAPInt nearest;

727

728

729

730

731

733

734

735 unsigned k = 1;

737 for (unsigned j = k - 1; j < k; j--) {

738

741 nearest = round(mu);

742

744 gsOrth = gramSchmidt();

745 }

748

750 (delta - mu * mu) *

752

753 k += 1;

754 } else {

755

756

758 gsOrth = gramSchmidt();

759 k = k > 1 ? k - 1 : 1;

760 }

761 }

762 }

763

767 IntMatrix normalized(numRows, numColumns);

768

769 DynamicAPInt lcmDenoms = DynamicAPInt(1);

770 for (unsigned i = 0; i < numRows; i++) {

771

772 for (unsigned j = 0; j < numColumns; j++)

773 lcmDenoms = lcm(lcmDenoms, at(i, j).den);

774

775 for (unsigned j = 0; j < numColumns; j++)

776 normalized(i, j) = (at(i, j) * lcmDenoms).getAsInteger();

777 }

778 return normalized;

779 }

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.

Example usage: Print .12, 3.4, 56.7 preAlign = ".", minSpacing = 1, .12 .12 3.4 3....

Eliminates variable at the specified position using Fourier-Motzkin variable elimination.