Matrix.cpp
Go to the documentation of this file.
/* Frobby: Software for monomial ideal computations.
Copyright (C) 2010 University of Aarhus
Contact Bjarke Hammersholt Roune for license information (www.broune.com)
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program. If not, see http://www.gnu.org/licenses/.
*/
#include "stdinc.h"
#include "Matrix.h"
#include "BigIntVector.h"
#include "ColumnPrinter.h"
#include <utility>
#include <sstream>
namespace {
void makeColumnIntegralPrimitive(Matrix& mat, size_t col) {
ASSERT(col < mat.getColCount());
if (mat.getRowCount() == 0)
return;
// Obtain gcd of numerators and lcm of denominators.
mpz_class numGcd = mat(0, col).get_num();
mpz_class denLcm = mat(0, col).get_den();
for (size_t row = 1; row < mat.getRowCount(); ++row) {
mpz_gcd(numGcd.get_mpz_t(), numGcd.get_mpz_t(),
mat(row, col).get_num_mpz_t());
mpz_lcm(denLcm.get_mpz_t(), denLcm.get_mpz_t(),
mat(row, col).get_den_mpz_t());
}
ASSERT(numGcd > 0);
ASSERT(denLcm > 0);
for (size_t row = 0; row < mat.getRowCount(); ++row) {
mat(row, col) /= numGcd;
mat(row, col) *= denLcm;
ASSERT(mat(row, col).get_den() == 1);
}
}
}
Matrix::Matrix(size_t rowCount, size_t colCount):
_rowCount(rowCount), _colCount(colCount), _entries(rowCount * colCount) {
}
void Matrix::resize(size_t rowCount, size_t colCount) {
if (rowCount == getRowCount() && colCount == getColCount())
return;
Matrix mat(rowCount, colCount);
size_t copyRowCount = std::min(rowCount, getRowCount());
size_t copyColCount = std::min(colCount, getColCount());
for (size_t row = 0; row < copyRowCount; ++row)
for (size_t col = 0; col < copyColCount; ++col)
mat(row, col) = (*this)(row, col);
swap(mat);
}
void Matrix::swap(Matrix& mat) {
using std::swap;
_entries.swap(mat._entries);
swap(_rowCount, mat._rowCount);
swap(_colCount, mat._colCount);
}
bool operator==(const Matrix& a, const Matrix& b) {
if (a.getRowCount() != b.getRowCount() ||
a.getColCount() != b.getColCount())
return false;
for (size_t row = 0; row < a.getRowCount(); ++row)
for (size_t col = 0; col < a.getColCount(); ++col)
if (a(row, col) != b(row, col))
return false;
return true;
}
ostream& operator<<(ostream& out, const Matrix& mat) {
ColumnPrinter pr;
print(pr, mat);
pr.print(out);
return out;
}
void print(FILE* file, const Matrix& mat) {
ostringstream out;
out << mat;
fputs(out.str().c_str(), file);
}
void print(ColumnPrinter& pr, const Matrix& mat) {
size_t baseCol = pr.getColumnCount();
for (size_t i = 0; i < mat.getColCount(); ++i)
pr.addColumn(false);
for (size_t col = 0; col < mat.getColCount(); ++col)
for (size_t row = 0; row < mat.getRowCount(); ++row)
pr[baseCol + col] << mat(row, col) << '\n';
}
void product(Matrix& prod, const Matrix& a, const Matrix& b) {
ASSERT(a.getColCount() == b.getRowCount());
prod.resize(a.getRowCount(), b.getColCount());
for (size_t r = 0; r < a.getRowCount(); ++r) {
for (size_t c = 0; c < b.getColCount(); ++c) {
prod(r, c) = 0;
for (size_t i = 0; i < a.getColCount(); ++i)
prod(r, c) += a(r, i) * b(i, c);
}
}
}
void transpose(Matrix& trans, const Matrix& mat) {
if (&trans == &mat) {
transpose(trans);
return;
}
trans.resize(mat.getColCount(), mat.getRowCount());
for (size_t row = 0; row < mat.getRowCount(); ++row)
for (size_t col = 0; col < mat.getColCount(); ++col)
trans(col, row) = mat(row, col);
}
void transpose(Matrix& mat) {
Matrix tmp(mat);
transpose(mat, tmp);
}
void addMultiplyRow(Matrix& mat, size_t resultRow,
size_t sourceRow, const mpq_class& mult) {
ASSERT(resultRow < mat.getRowCount());
ASSERT(sourceRow < mat.getRowCount());
for (size_t col = 0; col < mat.getColCount(); ++col)
mat(resultRow, col) += mat(sourceRow, col) * mult;
}
void multiplyRow(Matrix& mat, size_t row, const mpq_class& mult) {
for (size_t col = 0; col < mat.getColCount(); ++col)
mat(row, col) *= mult;
}
void swapRows(Matrix& mat, size_t row1, size_t row2) {
ASSERT(row1 < mat.getRowCount());
ASSERT(row2 < mat.getRowCount());
for (size_t col = 0; col < mat.getColCount(); ++col)
swap(mat(row1, col), mat(row2, col));
}
bool rowReduce(Matrix& mat) {
bool permutationOdd = false;
size_t rowsDone = 0;
for (size_t pivotCol = 0; pivotCol < mat.getColCount(); ++pivotCol) {
size_t pivotRow = rowsDone;
for (; pivotRow < mat.getRowCount(); ++pivotRow)
if (mat(pivotRow, pivotCol) != 0)
break;
if (pivotRow == mat.getRowCount())
continue;
if (rowsDone != pivotRow) {
permutationOdd = !permutationOdd;
swapRows(mat, rowsDone, pivotRow);
}
pivotRow = rowsDone;
++rowsDone;
for (size_t row = pivotRow + 1; row < mat.getRowCount(); ++row) {
if (row != pivotRow && mat(row, pivotCol) != 0) {
addMultiplyRow(mat, row, pivotRow,
-mat(row, pivotCol) / mat(pivotRow, pivotCol));
ASSERT(mat(row, pivotCol) == 0);
}
}
}
return permutationOdd;
}
void rowReduceFully(Matrix& mat) {
rowReduce(mat);
size_t pivotCol = 0;
size_t pivotRow = 0;
while (pivotRow < mat.getRowCount() &&
pivotCol < mat.getColCount()) {
if (mat(pivotRow, pivotCol) == 0) {
++pivotCol;
} else {
multiplyRow(mat, pivotRow, 1 / mat(pivotRow, pivotCol));
ASSERT(mat(pivotRow, pivotCol) == 1);
for (size_t row = 0; row < pivotRow; ++row) {
if (row != pivotRow && mat(row, pivotCol) != 0) {
addMultiplyRow(mat, row, pivotRow, -mat(row, pivotCol));
ASSERT(mat(row, pivotCol) == 0);
}
}
++pivotRow;
}
}
}
void subMatrix(Matrix& sub, const Matrix& mat,
size_t rowBegin, size_t rowEnd,
size_t colBegin, size_t colEnd) {
ASSERT(rowBegin <= rowEnd);
ASSERT(rowEnd <= mat.getRowCount());
ASSERT(colBegin <= colEnd);
ASSERT(colEnd <= mat.getColCount());
if (&sub == &mat) {
Matrix tmp;
subMatrix(tmp, mat, rowBegin, rowEnd, colBegin, colEnd);
sub.swap(tmp);
return;
}
sub.resize(rowEnd - rowBegin, colEnd - colBegin);
for (size_t row = rowBegin; row < rowEnd; ++row)
for (size_t col = colBegin; col < colEnd; ++col)
sub(row - rowBegin, col - colBegin) = mat(row, col);
}
void copyRow(Matrix& target, size_t targetRow,
const Matrix& source, size_t sourceRow) {
ASSERT(target.getColCount() == source.getColCount());
ASSERT(targetRow < target.getRowCount());
ASSERT(sourceRow < source.getRowCount());
size_t colCount = target.getColCount();
for (size_t col = 0; col < colCount; ++col)
target(targetRow, col) = source(sourceRow, col);
}
bool inverse(Matrix& inv, const Matrix& mat) {
ASSERT(mat.getRowCount() == mat.getColCount());
size_t size = mat.getRowCount();
inv = mat;
// Append identity matrix
inv.resize(size, size + size);
for (size_t i = 0; i < size; ++i)
inv(i, size + i) = 1;
rowReduceFully(inv);
if (inv(size - 1, size - 1) == 0)
return false; // not invertible
subMatrix(inv, inv, 0, size, size, 2 * size);
return true;
}
size_t matrixRank(const Matrix& matParam) {
Matrix mat(matParam);
rowReduceFully(mat);
// Find pivots
size_t rank = 0;
size_t col = 0;
size_t row = 0;
while (row < mat.getRowCount() && col < mat.getColCount()) {
if (mat(row, col) == 0) {
++col;
} else {
++rank;
++row;
}
}
return rank;
}
void nullSpace(Matrix& basis, const Matrix& matParam) {
Matrix mat(matParam);
rowReduceFully(mat);
// Find pivots
size_t rank = 0;
vector<char> colHasPivot(mat.getColCount());
{
size_t col = 0;
size_t row = 0;
while (row < mat.getRowCount() && col < mat.getColCount()) {
if (mat(row, col) == 0) {
++col;
} else {
++rank;
colHasPivot[col] = true;
++row;
}
}
}
// Construct basis
basis.resize(mat.getColCount(), mat.getColCount() - rank);
size_t nullCol = 0;
for (size_t col = 0; col < mat.getColCount(); ++col) {
ASSERT(nullCol <= basis.getColCount());
if (colHasPivot[col])
continue;
ASSERT(nullCol < basis.getColCount());
size_t row = 0;
for (size_t nullRow = 0; nullRow < basis.getRowCount(); ++nullRow) {
if (colHasPivot[nullRow]) {
basis(nullRow, nullCol) = -mat(row, col);
++row;
} else if (nullRow == col)
basis(nullRow, nullCol) = 1;
else
basis(nullRow, nullCol) = 0;
}
++nullCol;
}
ASSERT(nullCol == basis.getColCount());
// Make basis integer
for (size_t col = 0; col < basis.getColCount(); ++col)
makeColumnIntegralPrimitive(basis, col);
}
bool solve(Matrix& sol, const Matrix& lhs, const Matrix& rhs) {
ASSERT(lhs.getRowCount() == rhs.getRowCount());
// Append lhs|rhs and reduce
Matrix system = lhs;
system.resize(system.getRowCount(), system.getColCount() + rhs.getColCount());
size_t midCol = lhs.getColCount();
for (size_t col = 0; col < rhs.getColCount(); ++col)
for (size_t row = 0; row < rhs.getRowCount(); ++row)
system(row, midCol + col) = rhs(row, col);
rowReduceFully(system);
// Check if system has a solution
for (size_t row = 0; row < system.getRowCount(); ++row) {
for (size_t col = 0; col < midCol; ++col)
if (system(row, col) != 0)
goto hasLeftPivot;
for (size_t col = midCol; col < system.getColCount(); ++col)
if (system(row, col) != 0)
return false;
hasLeftPivot:;
}
// Extract solution
sol.resize(lhs.getColCount(), rhs.getColCount());
size_t row = 0;
for (size_t col = 0; col < midCol; ++col) {
if (row == system.getRowCount() || system(row, col) == 0) {
for (size_t r = 0; r < sol.getColCount(); ++r)
sol(col, r) = 0;
} else {
ASSERT(system(row, col) == 1);
for (size_t r = 0; r < sol.getColCount(); ++r)
sol(col, r) = system(row, midCol + r);
++row;
}
}
return true;
}
bool hasSameRowSpace(const Matrix& a, const Matrix& b) {
Matrix trA;
transpose(trA, a);
Matrix trB;
transpose(trB, b);
return hasSameColSpace(trA, trB);
}
bool hasSameColSpace(const Matrix& a, const Matrix& b) {
if (a.getRowCount() != b.getRowCount())
return false;
// A single row reduction of each of a and b would be a little more
// efficient.
Matrix dummy;
return solve(dummy, a, b) && solve(dummy, b, a);
}
mpq_class determinant(const Matrix& mat) {
ASSERT(mat.getRowCount() == mat.getColCount());
Matrix reduced(mat);
bool permutationOdd = rowReduce(reduced);
mpq_class det = permutationOdd ? -1 : 1;
for (size_t i = 0; i < reduced.getRowCount(); ++i)
det *= reduced(i, i);
return det;
}
namespace {
size_t getOppositeZeroRow(const Matrix& mat) {
// Let a,d and b,c be opposite vertices in a parallelogram. Then
// and only then b + c == d + a. We return the index of the row opposite
// to row 0. If the rows of mat are not the vertices of a parallelogram
// then we return mat.getRowCount().
if (mat.getRowCount() != 4)
return mat.getRowCount();
mpq_class tmp;
for (size_t opposite = 1; opposite < 4; ++opposite) {
bool isPara = true;
for (size_t col = 0; col < mat.getColCount(); ++col) {
tmp = mat(0, col) + mat(opposite, col);
for (size_t row = 1; row < 4; ++row)
if (row != opposite)
tmp -= mat(row, col);
if (tmp != 0) {
isPara = false;
break;
}
}
if (isPara)
return opposite;
}
return mat.getRowCount();
}
}
bool isParallelogram(const Matrix& mat) {
return getOppositeZeroRow(mat) != mat.getRowCount();
}
mpq_class getParallelogramAreaSq(const Matrix& mat) {
ASSERT(isParallelogram(mat));
size_t opposite = getOppositeZeroRow(mat);
size_t a;
for (a = 1; a < 4; ++a)
if (a != opposite)
break;
ASSERT(a < 4);
size_t b;
for (b = a + 1; b < 4; ++b)
if (b != opposite)
break;
ASSERT(b < 4);
// Translate to zero and drop the zero and sum vertices.
Matrix tmp(2, mat.getColCount());
for (size_t col = 0; col < mat.getColCount(); ++col) {
tmp(0, col) = mat(a, col) - mat(0, col);
tmp(1, col) = mat(b, col) - mat(0, col);
}
// Now the square of the area is det(tmp*transpose(tmp)).
Matrix trans;
transpose(trans, tmp);
Matrix prod;
product(prod, tmp, trans);
return determinant(prod);
}
BigIntVector.h
ColumnPrinter.h
swap
void swap(HilbertSlice &a, HilbertSlice &b)
Definition HilbertSlice.h:63
matrixRank
size_t matrixRank(const Matrix &matParam)
Returns the rank of mat.
Definition Matrix.cpp:276
multiplyRow
void multiplyRow(Matrix &mat, size_t row, const mpq_class &mult)
Multiplies row row with mult.
Definition Matrix.cpp:156
solve
bool solve(Matrix &sol, const Matrix &lhs, const Matrix &rhs)
Sets sol to some matrix such that lhs*sol=rhs and returns true if such a matrix exists.
Definition Matrix.cpp:348
getParallelogramAreaSq
mpq_class getParallelogramAreaSq(const Matrix &mat)
Returns the square of the area of the parallelogram whose vertices are the 4 rows of mat.
Definition Matrix.cpp:456
product
void product(Matrix &prod, const Matrix &a, const Matrix &b)
Sets prod to a * b.
Definition Matrix.cpp:116
isParallelogram
bool isParallelogram(const Matrix &mat)
Returns true if the rows of mat are the (4) vertices of a parallelogram.
Definition Matrix.cpp:452
addMultiplyRow
void addMultiplyRow(Matrix &mat, size_t resultRow, size_t sourceRow, const mpq_class &mult)
Adds mult times row sourceRow to row resultRow of mat.
Definition Matrix.cpp:147
rowReduce
bool rowReduce(Matrix &mat)
Reduces mat to row-echelon form, i.e.
Definition Matrix.cpp:169
swapRows
void swapRows(Matrix &mat, size_t row1, size_t row2)
Swaps row row1 and row row2 of mat.
Definition Matrix.cpp:161
nullSpace
void nullSpace(Matrix &basis, const Matrix &matParam)
Sets the columns of basis to a basis of the null space of mat.
Definition Matrix.cpp:296
determinant
mpq_class determinant(const Matrix &mat)
Returns the determinant of mat.
Definition Matrix.cpp:410
print
void print(FILE *file, const Matrix &mat)
Definition Matrix.cpp:101
subMatrix
void subMatrix(Matrix &sub, const Matrix &mat, size_t rowBegin, size_t rowEnd, size_t colBegin, size_t colEnd)
Sets sub to the sub-matrix of mat with rows in the interval [rowBegin, rowEnd) and columns in the int...
Definition Matrix.cpp:225
rowReduceFully
void rowReduceFully(Matrix &mat)
Reduces mat to reduced row-echelon form, i.e.
Definition Matrix.cpp:200
copyRow
void copyRow(Matrix &target, size_t targetRow, const Matrix &source, size_t sourceRow)
Copies row sourceRow from source to row targetRow of target.
Definition Matrix.cpp:246
hasSameRowSpace
bool hasSameRowSpace(const Matrix &a, const Matrix &b)
Returns true if a and b have the same row space.
Definition Matrix.cpp:390
inverse
bool inverse(Matrix &inv, const Matrix &mat)
Sets inv to the inverse of mat.
Definition Matrix.cpp:257
operator==
bool operator==(const Matrix &a, const Matrix &b)
Definition Matrix.cpp:82
operator<<
ostream & operator<<(ostream &out, const Matrix &mat)
Definition Matrix.cpp:94
hasSameColSpace
bool hasSameColSpace(const Matrix &a, const Matrix &b)
Returns true if a and b have the same column space.
Definition Matrix.cpp:400
transpose
void transpose(Matrix &trans, const Matrix &mat)
Sets trans to the transpose of mat.
Definition Matrix.cpp:129
Matrix.h
ColumnPrinter
Definition ColumnPrinter.h:27
ColumnPrinter::print
void print(ostream &out) const
Definition ColumnPrinter.cpp:70
ColumnPrinter::getColumnCount
size_t getColumnCount() const
Definition ColumnPrinter.cpp:61
ColumnPrinter::addColumn
void addColumn(bool flushLeft=true, const string &prefix=" ", const string &suffix="")
Definition ColumnPrinter.cpp:51
Matrix
Definition Matrix.h:26
Matrix::_colCount
size_t _colCount
Definition Matrix.h:57
Matrix::resize
void resize(size_t rowCount, size_t colCount)
Set the number of rows and columns.
Definition Matrix.cpp:61
Matrix::_entries
vector< mpq_class > _entries
Definition Matrix.h:58
Matrix::swap
void swap(Matrix &mat)
Definition Matrix.cpp:74
Matrix::getColCount
size_t getColCount() const
Definition Matrix.h:31
Matrix::getRowCount
size_t getRowCount() const
Definition Matrix.h:30
Matrix::_rowCount
size_t _rowCount
Definition Matrix.h:56
Matrix::Matrix
Matrix(size_t rowCount=0, size_t colCount=0)
Definition Matrix.cpp:57
stdinc.h
This header file includes common definitions and is included as the first line of code in every imple...
ASSERT
#define ASSERT(X)
Definition stdinc.h:86