Frobby 0.9.5
Matrix.cpp
Go to the documentation of this file.
1/* Frobby: Software for monomial ideal computations.
2 Copyright (C) 2010 University of Aarhus
3 Contact Bjarke Hammersholt Roune for license information (www.broune.com)
4
5 This program is free software; you can redistribute it and/or modify
6 it under the terms of the GNU General Public License as published by
7 the Free Software Foundation; either version 2 of the License, or
8 (at your option) any later version.
9
10 This program is distributed in the hope that it will be useful,
11 but WITHOUT ANY WARRANTY; without even the implied warranty of
12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 GNU General Public License for more details.
14
15 You should have received a copy of the GNU General Public License
16 along with this program. If not, see http://www.gnu.org/licenses/.
17*/
18#include "stdinc.h"
19#include "Matrix.h"
20
21#include "BigIntVector.h"
22#include "ColumnPrinter.h"
23
24#include <utility>
25#include <sstream>
26
27namespace {
31 void makeColumnIntegralPrimitive(Matrix& mat, size_t col) {
32 ASSERT(col < mat.getColCount());
33 if (mat.getRowCount() == 0)
34 return;
35
36 // Obtain gcd of numerators and lcm of denominators.
37 mpz_class numGcd = mat(0, col).get_num();
38 mpz_class denLcm = mat(0, col).get_den();
39 for (size_t row = 1; row < mat.getRowCount(); ++row) {
40 mpz_gcd(numGcd.get_mpz_t(), numGcd.get_mpz_t(),
41 mat(row, col).get_num_mpz_t());
42 mpz_lcm(denLcm.get_mpz_t(), denLcm.get_mpz_t(),
43 mat(row, col).get_den_mpz_t());
44 }
45 ASSERT(numGcd > 0);
46 ASSERT(denLcm > 0);
47
48 for (size_t row = 0; row < mat.getRowCount(); ++row) {
49 mat(row, col) /= numGcd;
50 mat(row, col) *= denLcm;
51
52 ASSERT(mat(row, col).get_den() == 1);
53 }
54 }
55}
56
57Matrix::Matrix(size_t rowCount, size_t colCount):
58 _rowCount(rowCount), _colCount(colCount), _entries(rowCount * colCount) {
59}
60
61void Matrix::resize(size_t rowCount, size_t colCount) {
62 if (rowCount == getRowCount() && colCount == getColCount())
63 return;
64 Matrix mat(rowCount, colCount);
65
66 size_t copyRowCount = std::min(rowCount, getRowCount());
67 size_t copyColCount = std::min(colCount, getColCount());
68 for (size_t row = 0; row < copyRowCount; ++row)
69 for (size_t col = 0; col < copyColCount; ++col)
70 mat(row, col) = (*this)(row, col);
71 swap(mat);
72}
73
74void Matrix::swap(Matrix& mat) {
75 using std::swap;
76
77 _entries.swap(mat._entries);
80}
81
82bool operator==(const Matrix& a, const Matrix& b) {
83 if (a.getRowCount() != b.getRowCount() ||
84 a.getColCount() != b.getColCount())
85 return false;
86
87 for (size_t row = 0; row < a.getRowCount(); ++row)
88 for (size_t col = 0; col < a.getColCount(); ++col)
89 if (a(row, col) != b(row, col))
90 return false;
91 return true;
92}
93
94ostream& operator<<(ostream& out, const Matrix& mat) {
96 print(pr, mat);
97 pr.print(out);
98 return out;
99}
100
101void print(FILE* file, const Matrix& mat) {
102 ostringstream out;
103 out << mat;
104 fputs(out.str().c_str(), file);
105}
106
107void print(ColumnPrinter& pr, const Matrix& mat) {
108 size_t baseCol = pr.getColumnCount();
109 for (size_t i = 0; i < mat.getColCount(); ++i)
110 pr.addColumn(false);
111 for (size_t col = 0; col < mat.getColCount(); ++col)
112 for (size_t row = 0; row < mat.getRowCount(); ++row)
113 pr[baseCol + col] << mat(row, col) << '\n';
114}
115
116void product(Matrix& prod, const Matrix& a, const Matrix& b) {
117 ASSERT(a.getColCount() == b.getRowCount());
118
119 prod.resize(a.getRowCount(), b.getColCount());
120 for (size_t r = 0; r < a.getRowCount(); ++r) {
121 for (size_t c = 0; c < b.getColCount(); ++c) {
122 prod(r, c) = 0;
123 for (size_t i = 0; i < a.getColCount(); ++i)
124 prod(r, c) += a(r, i) * b(i, c);
125 }
126 }
127}
128
129void transpose(Matrix& trans, const Matrix& mat) {
130 if (&trans == &mat) {
131 transpose(trans);
132 return;
133 }
134
135 trans.resize(mat.getColCount(), mat.getRowCount());
136 for (size_t row = 0; row < mat.getRowCount(); ++row)
137 for (size_t col = 0; col < mat.getColCount(); ++col)
138 trans(col, row) = mat(row, col);
139}
140
141void transpose(Matrix& mat) {
142 Matrix tmp(mat);
143 transpose(mat, tmp);
144}
145
146
147void addMultiplyRow(Matrix& mat, size_t resultRow,
148 size_t sourceRow, const mpq_class& mult) {
149 ASSERT(resultRow < mat.getRowCount());
150 ASSERT(sourceRow < mat.getRowCount());
151
152 for (size_t col = 0; col < mat.getColCount(); ++col)
153 mat(resultRow, col) += mat(sourceRow, col) * mult;
154}
155
156void multiplyRow(Matrix& mat, size_t row, const mpq_class& mult) {
157 for (size_t col = 0; col < mat.getColCount(); ++col)
158 mat(row, col) *= mult;
159}
160
161void swapRows(Matrix& mat, size_t row1, size_t row2) {
162 ASSERT(row1 < mat.getRowCount());
163 ASSERT(row2 < mat.getRowCount());
164
165 for (size_t col = 0; col < mat.getColCount(); ++col)
166 swap(mat(row1, col), mat(row2, col));
167}
168
169bool rowReduce(Matrix& mat) {
170 bool permutationOdd = false;
171
172 size_t rowsDone = 0;
173 for (size_t pivotCol = 0; pivotCol < mat.getColCount(); ++pivotCol) {
174 size_t pivotRow = rowsDone;
175 for (; pivotRow < mat.getRowCount(); ++pivotRow)
176 if (mat(pivotRow, pivotCol) != 0)
177 break;
178 if (pivotRow == mat.getRowCount())
179 continue;
180
181 if (rowsDone != pivotRow) {
182 permutationOdd = !permutationOdd;
183 swapRows(mat, rowsDone, pivotRow);
184 }
185 pivotRow = rowsDone;
186 ++rowsDone;
187
188 for (size_t row = pivotRow + 1; row < mat.getRowCount(); ++row) {
189 if (row != pivotRow && mat(row, pivotCol) != 0) {
190 addMultiplyRow(mat, row, pivotRow,
191 -mat(row, pivotCol) / mat(pivotRow, pivotCol));
192 ASSERT(mat(row, pivotCol) == 0);
193 }
194 }
195 }
196
197 return permutationOdd;
198}
199
201 rowReduce(mat);
202
204 size_t pivotCol = 0;
205 size_t pivotRow = 0;
206 while (pivotRow < mat.getRowCount() &&
207 pivotCol < mat.getColCount()) {
208 if (mat(pivotRow, pivotCol) == 0) {
209 ++pivotCol;
210 } else {
211 multiplyRow(mat, pivotRow, 1 / mat(pivotRow, pivotCol));
212 ASSERT(mat(pivotRow, pivotCol) == 1);
213 for (size_t row = 0; row < pivotRow; ++row) {
214 if (row != pivotRow && mat(row, pivotCol) != 0) {
215 addMultiplyRow(mat, row, pivotRow, -mat(row, pivotCol));
216 ASSERT(mat(row, pivotCol) == 0);
217 }
218 }
219
220 ++pivotRow;
221 }
222 }
223}
224
225void subMatrix(Matrix& sub, const Matrix& mat,
226 size_t rowBegin, size_t rowEnd,
227 size_t colBegin, size_t colEnd) {
228 ASSERT(rowBegin <= rowEnd);
229 ASSERT(rowEnd <= mat.getRowCount());
230 ASSERT(colBegin <= colEnd);
231 ASSERT(colEnd <= mat.getColCount());
232
233 if (&sub == &mat) {
234 Matrix tmp;
235 subMatrix(tmp, mat, rowBegin, rowEnd, colBegin, colEnd);
236 sub.swap(tmp);
237 return;
238 }
239
240 sub.resize(rowEnd - rowBegin, colEnd - colBegin);
241 for (size_t row = rowBegin; row < rowEnd; ++row)
242 for (size_t col = colBegin; col < colEnd; ++col)
243 sub(row - rowBegin, col - colBegin) = mat(row, col);
244}
245
246void copyRow(Matrix& target, size_t targetRow,
247 const Matrix& source, size_t sourceRow) {
248 ASSERT(target.getColCount() == source.getColCount());
249 ASSERT(targetRow < target.getRowCount());
250 ASSERT(sourceRow < source.getRowCount());
251
252 size_t colCount = target.getColCount();
253 for (size_t col = 0; col < colCount; ++col)
254 target(targetRow, col) = source(sourceRow, col);
255}
256
257bool inverse(Matrix& inv, const Matrix& mat) {
258 ASSERT(mat.getRowCount() == mat.getColCount());
259 size_t size = mat.getRowCount();
260
261 inv = mat;
262
263 // Append identity matrix
264 inv.resize(size, size + size);
265 for (size_t i = 0; i < size; ++i)
266 inv(i, size + i) = 1;
267
268 rowReduceFully(inv);
269 if (inv(size - 1, size - 1) == 0)
270 return false; // not invertible
271
272 subMatrix(inv, inv, 0, size, size, 2 * size);
273 return true;
274}
275
276size_t matrixRank(const Matrix& matParam) {
277 Matrix mat(matParam);
278 rowReduceFully(mat);
279
280 // Find pivots
281 size_t rank = 0;
282 size_t col = 0;
283 size_t row = 0;
284 while (row < mat.getRowCount() && col < mat.getColCount()) {
285 if (mat(row, col) == 0) {
286 ++col;
287 } else {
288 ++rank;
289 ++row;
290 }
291 }
292
293 return rank;
294}
295
296void nullSpace(Matrix& basis, const Matrix& matParam) {
297 Matrix mat(matParam);
298 rowReduceFully(mat);
299
300 // Find pivots
301 size_t rank = 0;
302 vector<char> colHasPivot(mat.getColCount());
303 {
304 size_t col = 0;
305 size_t row = 0;
306 while (row < mat.getRowCount() && col < mat.getColCount()) {
307 if (mat(row, col) == 0) {
308 ++col;
309 } else {
310 ++rank;
311 colHasPivot[col] = true;
312 ++row;
313 }
314 }
315 }
316
317 // Construct basis
318 basis.resize(mat.getColCount(), mat.getColCount() - rank);
319 size_t nullCol = 0;
320 for (size_t col = 0; col < mat.getColCount(); ++col) {
321 ASSERT(nullCol <= basis.getColCount());
322
323 if (colHasPivot[col])
324 continue;
325
326 ASSERT(nullCol < basis.getColCount());
327
328 size_t row = 0;
329 for (size_t nullRow = 0; nullRow < basis.getRowCount(); ++nullRow) {
330 if (colHasPivot[nullRow]) {
331 basis(nullRow, nullCol) = -mat(row, col);
332 ++row;
333 } else if (nullRow == col)
334 basis(nullRow, nullCol) = 1;
335 else
336 basis(nullRow, nullCol) = 0;
337 }
338
339 ++nullCol;
340 }
341 ASSERT(nullCol == basis.getColCount());
342
343 // Make basis integer
344 for (size_t col = 0; col < basis.getColCount(); ++col)
345 makeColumnIntegralPrimitive(basis, col);
346}
347
348bool solve(Matrix& sol, const Matrix& lhs, const Matrix& rhs) {
349 ASSERT(lhs.getRowCount() == rhs.getRowCount());
350
351 // Append lhs|rhs and reduce
352 Matrix system = lhs;
353 system.resize(system.getRowCount(), system.getColCount() + rhs.getColCount());
354 size_t midCol = lhs.getColCount();
355 for (size_t col = 0; col < rhs.getColCount(); ++col)
356 for (size_t row = 0; row < rhs.getRowCount(); ++row)
357 system(row, midCol + col) = rhs(row, col);
358
359 rowReduceFully(system);
360
361 // Check if system has a solution
362 for (size_t row = 0; row < system.getRowCount(); ++row) {
363 for (size_t col = 0; col < midCol; ++col)
364 if (system(row, col) != 0)
365 goto hasLeftPivot;
366 for (size_t col = midCol; col < system.getColCount(); ++col)
367 if (system(row, col) != 0)
368 return false;
369
370 hasLeftPivot:;
371 }
372
373 // Extract solution
374 sol.resize(lhs.getColCount(), rhs.getColCount());
375 size_t row = 0;
376 for (size_t col = 0; col < midCol; ++col) {
377 if (row == system.getRowCount() || system(row, col) == 0) {
378 for (size_t r = 0; r < sol.getColCount(); ++r)
379 sol(col, r) = 0;
380 } else {
381 ASSERT(system(row, col) == 1);
382 for (size_t r = 0; r < sol.getColCount(); ++r)
383 sol(col, r) = system(row, midCol + r);
384 ++row;
385 }
386 }
387 return true;
388}
389
390bool hasSameRowSpace(const Matrix& a, const Matrix& b) {
391 Matrix trA;
392 transpose(trA, a);
393
394 Matrix trB;
395 transpose(trB, b);
396
397 return hasSameColSpace(trA, trB);
398}
399
400bool hasSameColSpace(const Matrix& a, const Matrix& b) {
401 if (a.getRowCount() != b.getRowCount())
402 return false;
403
404 // A single row reduction of each of a and b would be a little more
405 // efficient.
406 Matrix dummy;
407 return solve(dummy, a, b) && solve(dummy, b, a);
408}
409
410mpq_class determinant(const Matrix& mat) {
411 ASSERT(mat.getRowCount() == mat.getColCount());
412
413 Matrix reduced(mat);
414 bool permutationOdd = rowReduce(reduced);
415
416 mpq_class det = permutationOdd ? -1 : 1;
417 for (size_t i = 0; i < reduced.getRowCount(); ++i)
418 det *= reduced(i, i);
419 return det;
420}
421
422namespace {
423 size_t getOppositeZeroRow(const Matrix& mat) {
424 // Let a,d and b,c be opposite vertices in a parallelogram. Then
425 // and only then b + c == d + a. We return the index of the row opposite
426 // to row 0. If the rows of mat are not the vertices of a parallelogram
427 // then we return mat.getRowCount().
428
429 if (mat.getRowCount() != 4)
430 return mat.getRowCount();
431
432 mpq_class tmp;
433 for (size_t opposite = 1; opposite < 4; ++opposite) {
434 bool isPara = true;
435 for (size_t col = 0; col < mat.getColCount(); ++col) {
436 tmp = mat(0, col) + mat(opposite, col);
437 for (size_t row = 1; row < 4; ++row)
438 if (row != opposite)
439 tmp -= mat(row, col);
440 if (tmp != 0) {
441 isPara = false;
442 break;
443 }
444 }
445 if (isPara)
446 return opposite;
447 }
448 return mat.getRowCount();
449 }
450}
451
452bool isParallelogram(const Matrix& mat) {
453 return getOppositeZeroRow(mat) != mat.getRowCount();
454}
455
456mpq_class getParallelogramAreaSq(const Matrix& mat) {
458 size_t opposite = getOppositeZeroRow(mat);
459
460 size_t a;
461 for (a = 1; a < 4; ++a)
462 if (a != opposite)
463 break;
464 ASSERT(a < 4);
465
466 size_t b;
467 for (b = a + 1; b < 4; ++b)
468 if (b != opposite)
469 break;
470 ASSERT(b < 4);
471
472 // Translate to zero and drop the zero and sum vertices.
473 Matrix tmp(2, mat.getColCount());
474 for (size_t col = 0; col < mat.getColCount(); ++col) {
475 tmp(0, col) = mat(a, col) - mat(0, col);
476 tmp(1, col) = mat(b, col) - mat(0, col);
477 }
478
479 // Now the square of the area is det(tmp*transpose(tmp)).
480 Matrix trans;
481 transpose(trans, tmp);
482 Matrix prod;
483 product(prod, tmp, trans);
484
485 return determinant(prod);
486}
void swap(HilbertSlice &a, HilbertSlice &b)
size_t matrixRank(const Matrix &matParam)
Returns the rank of mat.
Definition Matrix.cpp:276
void multiplyRow(Matrix &mat, size_t row, const mpq_class &mult)
Multiplies row row with mult.
Definition Matrix.cpp:156
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
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
void product(Matrix &prod, const Matrix &a, const Matrix &b)
Sets prod to a * b.
Definition Matrix.cpp:116
bool isParallelogram(const Matrix &mat)
Returns true if the rows of mat are the (4) vertices of a parallelogram.
Definition Matrix.cpp:452
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
bool rowReduce(Matrix &mat)
Reduces mat to row-echelon form, i.e.
Definition Matrix.cpp:169
void swapRows(Matrix &mat, size_t row1, size_t row2)
Swaps row row1 and row row2 of mat.
Definition Matrix.cpp:161
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
mpq_class determinant(const Matrix &mat)
Returns the determinant of mat.
Definition Matrix.cpp:410
void print(FILE *file, const Matrix &mat)
Definition Matrix.cpp:101
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
void rowReduceFully(Matrix &mat)
Reduces mat to reduced row-echelon form, i.e.
Definition Matrix.cpp:200
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
bool hasSameRowSpace(const Matrix &a, const Matrix &b)
Returns true if a and b have the same row space.
Definition Matrix.cpp:390
bool inverse(Matrix &inv, const Matrix &mat)
Sets inv to the inverse of mat.
Definition Matrix.cpp:257
bool operator==(const Matrix &a, const Matrix &b)
Definition Matrix.cpp:82
ostream & operator<<(ostream &out, const Matrix &mat)
Definition Matrix.cpp:94
bool hasSameColSpace(const Matrix &a, const Matrix &b)
Returns true if a and b have the same column space.
Definition Matrix.cpp:400
void transpose(Matrix &trans, const Matrix &mat)
Sets trans to the transpose of mat.
Definition Matrix.cpp:129
void print(ostream &out) const
size_t getColumnCount() const
void addColumn(bool flushLeft=true, const string &prefix=" ", const string &suffix="")
size_t _colCount
Definition Matrix.h:57
void resize(size_t rowCount, size_t colCount)
Set the number of rows and columns.
Definition Matrix.cpp:61
vector< mpq_class > _entries
Definition Matrix.h:58
void swap(Matrix &mat)
Definition Matrix.cpp:74
size_t getColCount() const
Definition Matrix.h:31
size_t getRowCount() const
Definition Matrix.h:30
size_t _rowCount
Definition Matrix.h:56
Matrix(size_t rowCount=0, size_t colCount=0)
Definition Matrix.cpp:57
This header file includes common definitions and is included as the first line of code in every imple...
#define ASSERT(X)
Definition stdinc.h:86