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));
49void 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);
82 double multiplier = +1;
83 for (
int row = 0; row < m_rows; row++) {
93int Matrix::fold2dIndexes (
int row,
int col)
const
95 return row * m_cols + col;
100 int foldedIndex = fold2dIndexes (row, col);
101 return m_vector [foldedIndex];
104void 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) {
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) {
139 return inverseGaussianElimination (significantDigits,
144 double epsilonThreshold)
const
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;
168 Matrix adjoint = cofactor.transpose ();
171 if (valueFailsEpsilonTest (determ,
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,
190 inv.set (0, 0, 1.0 / denominator);
196Matrix Matrix::inverseGaussianElimination (
int significantDigits,
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 (qAbs (working.get (rowFrom, colFirstWithNonZero)) <= 0) {
256 working.normalizeRow (rowFrom, colFirstWithNonZero, significantDigits, matrixConsistent);
262 for (rowTo = rowFrom + 1; rowTo <
rows (); rowTo++) {
264 if (qAbs (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;
280 working.normalizeRow (rowFrom, colFirstWithNonZero, significantDigits, matrixConsistentIter);
286 for (rowTo = rowFrom - 1; rowTo >= 0; rowTo--) {
288 if (qAbs (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);
314unsigned int Matrix::leadingZeros (
int row)
const
316 unsigned int sum = 0;
317 for (
int col = 0; col <
cols (); col++) {
318 if (qAbs (
get (row, col)) > 0) {
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));
352void Matrix::normalizeRow (
int rowToNormalize,
354 int significantDigits,
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,
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);
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);
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;
433void 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));
482bool Matrix::valueFailsEpsilonTest (
double value,
483 double epsilonThreshold)
const
485 return qAbs (value) < qAbs (epsilonThreshold);
#define ENGAUGE_ASSERT(cond)
Drop in replacement for Q_ASSERT.
MatrixConsistent
Indicates if matrix is consistent (i.e. has at least one solution)
Matrix class that supports arbitrary NxN size.
double determinant() const
Return the determinant of this matrix.
int rows() const
Height of 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....
void set(int row, int col, double value)
Set (row, col) element.
Matrix transpose() const
Return the transpose of the current matrix.
Matrix(int N)
Simple constructor of square matrix with initialization to identity matrix.
Matrix operator*(const Matrix &other) const
Multiplication operator with a matrix.
double get(int row, int col) const
Return (row, col) element.
QString toString() const
Dump matrix to a string.
int cols() const
Width of matrix.