7 #include "EngaugeAssert.h" 10 #include <QTextStream> 25 m_rows = other.
rows();
26 m_cols = other.
cols();
27 m_vector.resize (m_rows * m_cols);
28 for (
int row = 0; row < m_rows; row++) {
29 for (
int col = 0; col < m_cols; col++) {
30 set (row, col, other.
get (row, col));
37 m_rows = other.
rows();
38 m_cols = other.
cols();
39 m_vector.resize (m_rows * m_cols);
40 for (
int row = 0; row < m_rows; row++) {
41 for (
int col = 0; col < m_cols; col++) {
42 set (row, col, other.
get (row, col));
49 void Matrix::addRowToAnotherWithScaling (
int rowFrom,
53 for (
int col = 0; col <
cols (); col++) {
54 double oldValueFrom =
get (rowFrom, col);
55 double oldValueTo =
get (rowTo, col);
56 double newValueTo = oldValueFrom * factor + oldValueTo;
57 set (rowTo, col, newValueTo);
68 ENGAUGE_ASSERT (m_rows == m_cols);
82 double multiplier = +1;
83 for (
int row = 0; row < m_rows; row++) {
85 rtn += multiplier *
get (row, COL) * min.
determinant ();
93 int Matrix::fold2dIndexes (
int row,
int col)
const 95 return row * m_cols + col;
100 int foldedIndex = fold2dIndexes (row, col);
101 return m_vector [foldedIndex];
104 void Matrix::initialize (
int rows,
110 for (
int row = 0; row < m_rows; row++) {
111 for (
int col = 0; col < m_cols; col++) {
114 if (row == col && m_rows == m_cols) {
124 MatrixConsistent &matrixConsistent)
const 128 for (
int row = 0; row < m_rows; row++) {
129 for (
int col = 0; col < m_cols; col++) {
130 double value = qAbs (
get (row, col));
131 if (value > maxValue) {
138 matrixConsistent = MATRIX_CONSISTENT;
139 return inverseGaussianElimination (significantDigits,
143 Matrix Matrix::inverseCramersRule (MatrixConsistent & matrixConsistent,
144 double epsilonThreshold)
const 146 ENGAUGE_ASSERT (m_rows == m_cols);
154 double multiplierStartForRow = 1.0;
156 for (row = 0; row < m_rows; row++) {
157 double multiplier = multiplierStartForRow;
158 for (col = 0; col < m_cols; col++) {
162 cofactor.set (row, col, element);
164 multiplierStartForRow *= -1.0;
171 if (valueFailsEpsilonTest (determ,
173 matrixConsistent = MATRIX_INCONSISTENT;
178 for (row = 0; row < m_rows; row++) {
179 for (col = 0; col < m_cols; col++) {
180 inv.set (row, col, adjoint.
get (row, col) / determ);
184 double denominator =
get (0, 0);
185 if (valueFailsEpsilonTest (denominator,
187 matrixConsistent = MATRIX_INCONSISTENT;
190 inv.set (0, 0, 1.0 / denominator);
196 Matrix Matrix::inverseGaussianElimination (
int significantDigits,
197 MatrixConsistent &matrixConsistent)
const 201 int row, col, rowFrom, rowTo;
205 QVector<int> rowIndexes (
rows ());
206 for (row = 0; row <
rows (); row++) {
207 rowIndexes [row] = row;
213 for (row = 0; row <
rows (); row++) {
214 for (col = 0; col <
cols (); col++) {
215 double value =
get (row, col);
216 working.set (row, col, value);
221 for (
int row1 = 0; row1 <
rows (); row1++) {
222 for (
int row2 = row1 + 1; row2 <
rows (); row2++) {
223 if (working.leadingZeros (row1) > working.leadingZeros (row2)) {
224 working.switchRows (row1, row2);
227 int row1Index = rowIndexes [row1];
228 int row2Index = rowIndexes [row2];
229 rowIndexes [row1] = row2Index;
230 rowIndexes [row2] = row1Index;
236 for (row = 0; row <
rows (); row++) {
237 for (col =
cols (); col < 2 *
cols (); col++) {
238 int colIdentitySubmatrix = col -
cols ();
239 double value = (row == colIdentitySubmatrix) ? 1 : 0;
240 working.set (row, col, value);
245 for (rowFrom = 0; rowFrom <
rows (); rowFrom++) {
246 int colFirstWithNonZero = rowFrom;
250 if (working.get (rowFrom, colFirstWithNonZero) == 0) {
251 matrixConsistent = MATRIX_INCONSISTENT;
256 working.normalizeRow (rowFrom, colFirstWithNonZero, significantDigits, matrixConsistent);
257 if (matrixConsistent == MATRIX_INCONSISTENT) {
262 for (rowTo = rowFrom + 1; rowTo <
rows (); rowTo++) {
264 if (working.get (rowTo, colFirstWithNonZero) != 0) {
267 double denominator = working.get (rowFrom, colFirstWithNonZero);
268 double factor = -1.0 * working.get (rowTo, colFirstWithNonZero) / denominator;
269 working.addRowToAnotherWithScaling (rowFrom, rowTo, factor);
275 for (rowFrom =
rows () - 1; rowFrom >= 0; rowFrom--) {
276 int colFirstWithNonZero = rowFrom;
279 MatrixConsistent matrixConsistent;
280 working.normalizeRow (rowFrom, colFirstWithNonZero, significantDigits, matrixConsistent);
281 if (matrixConsistent == MATRIX_INCONSISTENT) {
286 for (rowTo = rowFrom - 1; rowTo >= 0; rowTo--) {
288 if (working.get (rowTo, colFirstWithNonZero) != 0) {
291 double denominator = working.get (rowFrom, colFirstWithNonZero);
292 double factor = -1.0 * working.get (rowTo, colFirstWithNonZero) / denominator;
293 working.addRowToAnotherWithScaling (rowFrom, rowTo, factor);
300 for (row = 0; row < working.rows (); row++) {
302 int rowBeforeSort = rowIndexes [row];
304 for (col = 0; col < working.rows (); col++) {
305 int colRightHalf = col + inv.cols ();
306 double value = working.get (rowBeforeSort, colRightHalf);
307 inv.set (row, col, value);
314 unsigned int Matrix::leadingZeros (
int row)
const 316 unsigned int sum = 0;
317 for (
int col = 0; col <
cols (); col++) {
318 if (
get (row, col) != 0) {
328 ENGAUGE_ASSERT (m_rows == m_cols);
330 Matrix outMinor (m_rows - 1);
332 for (
int row = 0; row < m_rows; row++) {
334 if (row != rowOmit) {
337 for (
int col = 0; col < m_cols; col++) {
339 if (col != colOmit) {
341 outMinor.
set (rowMinor, colMinor,
get (row, col));
352 void Matrix::normalizeRow (
int rowToNormalize,
354 int significantDigits,
355 MatrixConsistent &matrixConsistent)
357 double denominator =
get (rowToNormalize, colToNormalize);
360 double smallestAbsValue = 0;
361 for (
int col = 0; col <
cols (); col++) {
362 double absValue = qAbs (
get (rowToNormalize, 0));
363 if (col == 0 || absValue < smallestAbsValue) {
364 smallestAbsValue = absValue;
367 double epsilonThreshold = smallestAbsValue / qPow (10.0, significantDigits);
369 if (valueFailsEpsilonTest (denominator,
372 matrixConsistent = MATRIX_INCONSISTENT;
376 matrixConsistent = MATRIX_CONSISTENT;
378 double factor = 1.0 / denominator;
379 for (
int col = 0; col <
cols (); col++) {
380 double value =
get (rowToNormalize, col);
381 set (rowToNormalize, col, factor * value);
388 ENGAUGE_ASSERT (m_cols == other.
rows ());
392 for (
int row = 0; row < m_rows; row++) {
393 for (
int col = 0; col < other.
cols (); col++) {
395 for (
int index = 0; index < m_cols; index++) {
396 sum +=
get (row, index) * other.
get (index, col);
398 out.set (row, col, sum);
407 ENGAUGE_ASSERT (m_cols == other.size ());
411 for (
int row = 0; row < m_rows; row++) {
413 for (
int col = 0; col < m_cols; col++) {
414 sum +=
get (row, col) * other [col];
430 m_vector [fold2dIndexes (row, col)] = value;
433 void Matrix::switchRows (
int row1,
437 for (
int col = 0; col <
cols (); col++) {
438 double temp1 =
get (row1, col);
439 double temp2 =
get (row2, col);
440 set (row2, col, temp2);
441 set (row1, col, temp1);
448 QTextStream str (&out);
451 for (
int row = 0; row <
rows (); row++) {
456 for (
int col = 0; col <
cols (); col++) {
460 str <<
get (row, col);
471 Matrix out (m_cols, m_rows);
473 for (
int row = 0; row < m_rows; row++) {
474 for (
int col = 0; col < m_cols; col++) {
475 out.
set (col, row,
get (row, col));
482 bool Matrix::valueFailsEpsilonTest (
double value,
483 double epsilonThreshold)
const 485 return qAbs (value) < qAbs (epsilonThreshold);
Matrix operator*(const Matrix &other) const
Multiplication operator with a matrix.
Matrix transpose() const
Return the transpose of the current matrix.
Matrix & operator=(const Matrix &matrix)
Assignment operator.
Matrix inverse(int significantDigits, MatrixConsistent &matrixConsistent) const
Return the inverse of this matrix.
Matrix minorReduced(int rowOmit, int colOmit) const
Return minor matrix which is the original with the specified row and column omitted. The name 'minor' is a reserved word.
double determinant() const
Return the determinant of this matrix.
Matrix(int N)
Simple constructor of square matrix with initialization to identity matrix.
int cols() const
Width of matrix.
Matrix class that supports arbitrary NxN size.
int rows() const
Height of matrix.
double get(int row, int col) const
Return (row, col) element.
void set(int row, int col, double value)
Set (row, col) element.
QString toString() const
Dump matrix to a string.