7 #include "EngaugeAssert.h"
10 #include <QTextStream>
20 initialize (rows, cols);
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 return m_vector [fold2dIndexes (row, col)];
103 void Matrix::initialize (
int rows,
108 m_vector.resize (rows * cols);
109 for (
int row = 0; row < m_rows; row++) {
110 for (
int col = 0; col < m_cols; col++) {
113 if (row == col && m_rows == m_cols) {
123 MatrixConsistent &matrixConsistent)
const
127 for (
int row = 0; row < m_rows; row++) {
128 for (
int col = 0; col < m_cols; col++) {
129 double value = qAbs (
get (row, col));
130 if (value > maxValue) {
136 double epsilonThreshold = maxValue / qPow (10.0, significantDigits);
139 matrixConsistent = MATRIX_CONSISTENT;
140 return inverseGaussianElimination (matrixConsistent,
144 Matrix Matrix::inverseCramersRule (MatrixConsistent & matrixConsistent,
145 double epsilonThreshold)
const
147 ENGAUGE_ASSERT (m_rows == m_cols);
155 double multiplierStartForRow = 1.0;
157 for (row = 0; row < m_rows; row++) {
158 double multiplier = multiplierStartForRow;
159 for (col = 0; col < m_cols; col++) {
163 cofactor.set (row, col, element);
165 multiplierStartForRow *= -1.0;
172 if (valueFailsEpsilonTest (determ,
174 matrixConsistent = MATRIX_INCONSISTENT;
179 for (row = 0; row < m_rows; row++) {
180 for (col = 0; col < m_cols; col++) {
181 inv.set (row, col, adjoint.
get (row, col) / determ);
185 double denominator =
get (0, 0);
186 if (valueFailsEpsilonTest (denominator,
188 matrixConsistent = MATRIX_INCONSISTENT;
191 inv.set (0, 0, 1.0 / denominator);
197 Matrix Matrix::inverseGaussianElimination (MatrixConsistent &matrixConsistent,
198 double epsilonThreshold)
const
202 int row, col, rowFrom, rowTo;
206 QVector<int> rowIndexes (rows ());
207 for (row = 0; row <
rows (); row++) {
208 rowIndexes [row] = row;
213 Matrix working (rows (), 2 * cols ());
214 for (row = 0; row <
rows (); row++) {
215 for (col = 0; col <
cols (); col++) {
216 double value =
get (row, col);
217 working.set (row, col, value);
222 for (
int row1 = 0; row1 <
rows (); row1++) {
223 for (
int row2 = row1 + 1; row2 <
rows (); row2++) {
224 if (working.leadingZeros (row1) > working.leadingZeros (row2)) {
225 working.switchRows (row1, row2);
228 int row1Index = rowIndexes [row1];
229 int row2Index = rowIndexes [row2];
230 rowIndexes [row1] = row2Index;
231 rowIndexes [row2] = row1Index;
237 for (row = 0; row <
rows (); row++) {
238 for (col = cols (); col < 2 *
cols (); col++) {
239 int colIdentitySubmatrix = col -
cols ();
240 double value = (row == colIdentitySubmatrix) ? 1 : 0;
241 working.set (row, col, value);
246 for (rowFrom = 0; rowFrom <
rows (); rowFrom++) {
247 int colFirstWithNonZero = rowFrom;
251 if (working.get (rowFrom, colFirstWithNonZero) == 0) {
252 matrixConsistent = MATRIX_INCONSISTENT;
257 working.normalizeRow (rowFrom, colFirstWithNonZero, matrixConsistent, epsilonThreshold);
258 if (matrixConsistent == MATRIX_INCONSISTENT) {
263 for (rowTo = rowFrom + 1; rowTo <
rows (); rowTo++) {
265 if (working.get (rowTo, colFirstWithNonZero) != 0) {
268 double denominator = working.get (rowFrom, colFirstWithNonZero);
269 if (valueFailsEpsilonTest (denominator, epsilonThreshold)) {
270 matrixConsistent = MATRIX_INCONSISTENT;
273 double factor = -1.0 * working.get (rowTo, colFirstWithNonZero) / denominator;
274 working.addRowToAnotherWithScaling (rowFrom, rowTo, factor);
280 for (rowFrom = rows () - 1; rowFrom >= 0; rowFrom--) {
281 int colFirstWithNonZero = rowFrom;
284 MatrixConsistent matrixConsistent;
285 working.normalizeRow (rowFrom, colFirstWithNonZero, matrixConsistent, epsilonThreshold);
286 if (matrixConsistent == MATRIX_INCONSISTENT) {
291 for (rowTo = rowFrom - 1; rowTo >= 0; rowTo--) {
293 if (working.get (rowTo, colFirstWithNonZero) != 0) {
296 double denominator = working.get (rowFrom, colFirstWithNonZero);
297 if (valueFailsEpsilonTest (denominator, epsilonThreshold)) {
298 matrixConsistent = MATRIX_INCONSISTENT;
301 double factor = -1.0 * working.get (rowTo, colFirstWithNonZero) / denominator;
302 working.addRowToAnotherWithScaling (rowFrom, rowTo, factor);
309 for (row = 0; row < working.rows (); row++) {
311 int rowBeforeSort = rowIndexes [row];
313 for (col = 0; col < working.rows (); col++) {
314 int colRightHalf = col + inv.cols ();
315 double value = working.get (rowBeforeSort, colRightHalf);
316 inv.set (row, col, value);
323 unsigned int Matrix::leadingZeros (
int row)
const
325 unsigned int sum = 0;
326 for (
int col = 0; col <
cols (); col++) {
327 if (
get (row, col) != 0) {
337 ENGAUGE_ASSERT (m_rows == m_cols);
339 Matrix outMinor (m_rows - 1);
341 for (
int row = 0; row < m_rows; row++) {
343 if (row != rowOmit) {
346 for (
int col = 0; col < m_cols; col++) {
348 if (col != colOmit) {
350 outMinor.
set (rowMinor, colMinor,
get (row, col));
361 void Matrix::normalizeRow (
int rowToNormalize,
363 MatrixConsistent &matrixConsistent,
364 double epsilonThreshold)
366 double denominator =
get (rowToNormalize, colToNormalize);
368 if (valueFailsEpsilonTest (denominator,
371 matrixConsistent = MATRIX_INCONSISTENT;
375 matrixConsistent = MATRIX_CONSISTENT;
377 double factor = 1.0 / denominator;
378 for (
int col = 0; col <
cols (); col++) {
379 double value =
get (rowToNormalize, col);
380 set (rowToNormalize, col, factor * value);
387 ENGAUGE_ASSERT (m_cols == other.
rows ());
391 for (
int row = 0; row < m_rows; row++) {
392 for (
int col = 0; col < other.
cols (); col++) {
394 for (
int index = 0; index < m_cols; index++) {
395 sum +=
get (row, index) * other.
get (index, col);
397 out.set (row, col, sum);
406 ENGAUGE_ASSERT (m_cols == other.size ());
410 for (
int row = 0; row < m_rows; row++) {
412 for (
int col = 0; col < m_cols; col++) {
413 sum +=
get (row, col) * other [col];
429 m_vector [fold2dIndexes (row, col)] = value;
432 void Matrix::switchRows (
int row1,
436 for (
int col = 0; col <
cols (); col++) {
437 double temp1 =
get (row1, col);
438 double temp2 =
get (row2, col);
439 set (row2, col, temp2);
440 set (row1, col, temp1);
447 QTextStream str (&out);
450 for (
int row = 0; row <
rows (); row++) {
455 for (
int col = 0; col <
cols (); col++) {
459 str <<
get (row, col);
470 Matrix out (m_cols, m_rows);
472 for (
int row = 0; row < m_rows; row++) {
473 for (
int col = 0; col < m_cols; col++) {
474 out.
set (col, row,
get (row, col));
481 bool Matrix::valueFailsEpsilonTest (
double value,
482 double epsilonThreshold)
const
484 return (qAbs (value) < qAbs (epsilonThreshold));
int rows() const
Height of 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.
Matrix inverse(int significantDigits, MatrixConsistent &matrixConsistent) const
Return the inverse of this matrix.
Matrix operator*(const Matrix &other) const
Multiplication operator with a matrix.
double determinant() const
Return the determinant of this matrix.
Matrix & operator=(const Matrix &matrix)
Assignment operator.
Matrix(int N)
Simple constructor of square matrix with initialization to identity matrix.
Matrix transpose() const
Return the transpose of the current matrix.
Matrix class that supports arbitrary NxN size.
QString toString() const
Dump matrix to a string.
int cols() const
Width of matrix.
void set(int row, int col, double value)
Set (row, col) element.
double get(int row, int col) const
Return (row, col) element.