7 #include "EngaugeAssert.h" 19 initialize (rows, cols);
24 m_rows = other.
rows();
25 m_cols = other.
cols();
26 m_vector.resize (m_rows * m_cols);
27 for (
int row = 0; row < m_rows; row++) {
28 for (
int col = 0; col < m_cols; col++) {
29 set (row, col, other.
get (row, col));
36 m_rows = other.
rows();
37 m_cols = other.
cols();
38 m_vector.resize (m_rows * m_cols);
39 for (
int row = 0; row < m_rows; row++) {
40 for (
int col = 0; col < m_cols; col++) {
41 set (row, col, other.
get (row, col));
48 void Matrix::addRowToAnotherWithScaling (
int rowFrom,
52 for (
int col = 0; col <
cols (); col++) {
53 double oldValueFrom =
get (rowFrom, col);
54 double oldValueTo =
get (rowTo, col);
55 double newValueTo = oldValueFrom * factor + oldValueTo;
56 set (rowTo, col, newValueTo);
67 ENGAUGE_ASSERT (m_rows == m_cols);
81 double multiplier = +1;
82 for (
int row = 0; row < m_rows; row++) {
84 rtn += multiplier *
get (row, COL) * min.
determinant ();
92 int Matrix::fold2dIndexes (
int row,
int col)
const 94 return row * m_cols + col;
99 return m_vector [fold2dIndexes (row, col)];
102 void Matrix::initialize (
int rows,
107 m_vector.resize (
rows * cols);
108 for (
int row = 0; row < m_rows; row++) {
109 for (
int col = 0; col < m_cols; col++) {
112 if (row == col && m_rows == m_cols) {
124 return inverseGaussianElimination ();
127 Matrix Matrix::inverseCramersRule ()
const 129 ENGAUGE_ASSERT (m_rows == m_cols);
137 double multiplierStartForRow = 1.0;
139 for (row = 0; row < m_rows; row++) {
140 double multiplier = multiplierStartForRow;
141 for (col = 0; col < m_cols; col++) {
145 cofactor.
set (row, col, element);
147 multiplierStartForRow *= -1.0;
156 for (row = 0; row < m_rows; row++) {
157 for (col = 0; col < m_cols; col++) {
158 inv.
set (row, col, adjoint.
get (row, col) / determ);
162 inv.
set (0, 0, 1.0 /
get (0, 0));
168 Matrix Matrix::inverseGaussianElimination ()
const 172 int row, col, rowFrom, rowTo;
175 QVector<int> rowIndexes (
rows ());
176 for (row = 0; row <
rows (); row++) {
177 rowIndexes [row] = row;
183 for (row = 0; row <
rows (); row++) {
184 for (col = 0; col <
cols (); col++) {
185 double value =
get (row, col);
186 working.
set (row, col, value);
191 for (
int row1 = 0; row1 <
rows (); row1++) {
192 for (
int row2 = row1 + 1; row2 <
rows (); row2++) {
193 if (working.leadingZeros (row1) > working.leadingZeros (row2)) {
194 working.switchRows (row1, row2);
197 int row1Index = rowIndexes [row1];
198 int row2Index = rowIndexes [row2];
199 rowIndexes [row1] = row2Index;
200 rowIndexes [row2] = row1Index;
206 for (row = 0; row <
rows (); row++) {
207 for (col =
cols (); col < 2 *
cols (); col++) {
208 int colIdentitySubmatrix = col -
cols ();
209 double value = (row == colIdentitySubmatrix) ? 1 : 0;
210 working.
set (row, col, value);
215 for (rowFrom = 0; rowFrom <
rows (); rowFrom++) {
216 int colFirstWithNonZero = rowFrom;
220 if (working.
get (rowFrom, colFirstWithNonZero) == 0) {
227 working.normalizeRow (rowFrom, colFirstWithNonZero);
230 for (rowTo = rowFrom + 1; rowTo <
rows (); rowTo++) {
232 if (working.
get (rowTo, colFirstWithNonZero) != 0) {
235 double factor = -1.0 * working.
get (rowTo, colFirstWithNonZero) / working.
get (rowFrom, colFirstWithNonZero);
236 working.addRowToAnotherWithScaling (rowFrom, rowTo, factor);
242 for (rowFrom = rows () - 1; rowFrom >= 0; rowFrom--) {
243 int colFirstWithNonZero = rowFrom;
246 working.normalizeRow (rowFrom, colFirstWithNonZero);
249 for (rowTo = rowFrom - 1; rowTo >= 0; rowTo--) {
251 if (working.
get (rowTo, colFirstWithNonZero) != 0) {
254 double factor = -1.0 * working.
get (rowTo, colFirstWithNonZero) / working.
get (rowFrom, colFirstWithNonZero);
255 working.addRowToAnotherWithScaling (rowFrom, rowTo, factor);
263 for (row = 0; row < working.
rows (); row++) {
265 int rowBeforeSort = rowIndexes [row];
267 for (col = 0; col < working.
rows (); col++) {
268 int colRightHalf = col + inv.
cols ();
269 double value = working.
get (rowBeforeSort, colRightHalf);
270 inv.set (row, col, value);
277 unsigned int Matrix::leadingZeros (
int row)
const 279 unsigned int sum = 0;
280 for (
int col = 0; col <
cols (); col++) {
281 if (
get (row, col) != 0) {
291 ENGAUGE_ASSERT (m_rows == m_cols);
293 Matrix outMinor (m_rows - 1);
295 for (
int row = 0; row < m_rows; row++) {
297 if (row != rowOmit) {
300 for (
int col = 0; col < m_cols; col++) {
302 if (col != colOmit) {
304 outMinor.
set (rowMinor, colMinor,
get (row, col));
315 void Matrix::normalizeRow (
int rowToNormalize,
318 ENGAUGE_ASSERT (
get (rowToNormalize, colToNormalize) != 0);
320 double factor = 1.0 /
get (rowToNormalize, colToNormalize);
321 for (
int col = 0; col <
cols (); col++) {
322 double value =
get (rowToNormalize, col);
323 set (rowToNormalize, col, factor * value);
329 ENGAUGE_ASSERT (m_cols == other.
rows ());
333 for (
int row = 0; row < m_rows; row++) {
334 for (
int col = 0; col < other.
cols (); col++) {
336 for (
int index = 0; index < m_cols; index++) {
337 sum +=
get (row, index) * other.
get (index, col);
339 out.set (row, col, sum);
348 ENGAUGE_ASSERT (m_cols == other.size ());
352 for (
int row = 0; row < m_rows; row++) {
354 for (
int col = 0; col < m_cols; col++) {
355 sum +=
get (row, col) * other [col];
371 m_vector [fold2dIndexes (row, col)] = value;
374 void Matrix::switchRows (
int row1,
378 for (
int col = 0; col <
cols (); col++) {
379 double temp1 =
get (row1, col);
380 double temp2 =
get (row2, col);
381 set (row2, col, temp2);
382 set (row1, col, temp1);
389 QTextStream str (&out);
392 for (
int row = 0; row <
rows (); row++) {
397 for (
int col = 0; col <
cols (); col++) {
401 str <<
get (row, col);
412 Matrix out (m_cols, m_rows);
414 for (
int row = 0; row < m_rows; row++) {
415 for (
int col = 0; col < m_cols; col++) {
416 out.
set (col, row,
get (row, col));
Matrix inverse() const
Return the inverse of this matrix.
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 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.