00001
00002
00003
00004
00005
00006
00007 #include "EngaugeAssert.h"
00008 #include "Matrix.h"
00009 #include <QTextStream>
00010
00011 Matrix::Matrix (int N)
00012 {
00013 initialize (N, N);
00014 }
00015
00016 Matrix::Matrix (int rows,
00017 int cols)
00018 {
00019 initialize (rows, cols);
00020 }
00021
00022 Matrix::Matrix (const Matrix &other)
00023 {
00024 m_rows = other.rows();
00025 m_cols = other.cols();
00026 m_vector.resize (m_rows * m_cols);
00027 for (int row = 0; row < m_rows; row++) {
00028 for (int col = 0; col < m_cols; col++) {
00029 set (row, col, other.get (row, col));
00030 }
00031 }
00032 }
00033
00034 Matrix &Matrix::operator= (const Matrix &other)
00035 {
00036 m_rows = other.rows();
00037 m_cols = other.cols();
00038 m_vector.resize (m_rows * m_cols);
00039 for (int row = 0; row < m_rows; row++) {
00040 for (int col = 0; col < m_cols; col++) {
00041 set (row, col, other.get (row, col));
00042 }
00043 }
00044
00045 return *this;
00046 }
00047
00048 void Matrix::addRowToAnotherWithScaling (int rowFrom,
00049 int rowTo,
00050 double factor)
00051 {
00052 for (int col = 0; col < cols (); col++) {
00053 double oldValueFrom = get (rowFrom, col);
00054 double oldValueTo = get (rowTo, col);
00055 double newValueTo = oldValueFrom * factor + oldValueTo;
00056 set (rowTo, col, newValueTo);
00057 }
00058 }
00059
00060 int Matrix::cols () const
00061 {
00062 return m_cols;
00063 }
00064
00065 double Matrix::determinant () const
00066 {
00067 ENGAUGE_ASSERT (m_rows == m_cols);
00068
00069 double rtn;
00070
00071 if (m_rows == 1) {
00072
00073 rtn = m_vector [0];
00074
00075 } else {
00076
00077 const int COL = 0;
00078
00079
00080 rtn = 0.0;
00081 double multiplier = +1;
00082 for (int row = 0; row < m_rows; row++) {
00083 Matrix min = minorReduced (row, COL);
00084 rtn += multiplier * get (row, COL) * min.determinant ();
00085 multiplier *= -1.0;
00086 }
00087 }
00088
00089 return rtn;
00090 }
00091
00092 int Matrix::fold2dIndexes (int row, int col) const
00093 {
00094 return row * m_cols + col;
00095 }
00096
00097 double Matrix::get (int row, int col) const
00098 {
00099 return m_vector [fold2dIndexes (row, col)];
00100 }
00101
00102 void Matrix::initialize (int rows,
00103 int cols)
00104 {
00105 m_rows = rows;
00106 m_cols = cols;
00107 m_vector.resize (rows * cols);
00108 for (int row = 0; row < m_rows; row++) {
00109 for (int col = 0; col < m_cols; col++) {
00110
00111
00112 if (row == col && m_rows == m_cols) {
00113 set (row, col, 1.0);
00114 } else {
00115 set (row, col, 0.0);
00116 }
00117 }
00118 }
00119 }
00120
00121 Matrix Matrix::inverse () const
00122 {
00123
00124 return inverseGaussianElimination ();
00125 }
00126
00127 Matrix Matrix::inverseCramersRule () const
00128 {
00129 ENGAUGE_ASSERT (m_rows == m_cols);
00130
00131 Matrix inv (m_rows);
00132 int row, col;
00133
00134 if (m_rows > 1) {
00135
00136
00137 double multiplierStartForRow = 1.0;
00138 Matrix cofactor (m_rows);
00139 for (row = 0; row < m_rows; row++) {
00140 double multiplier = multiplierStartForRow;
00141 for (col = 0; col < m_cols; col++) {
00142 Matrix min = minorReduced (row, col);
00143 double element = multiplier * min.determinant ();
00144 multiplier *= -1.0;
00145 cofactor.set (row, col, element);
00146 }
00147 multiplierStartForRow *= -1.0;
00148 }
00149
00150
00151 Matrix adjoint = cofactor.transpose ();
00152
00153 double determ = determinant ();
00154
00155
00156 for (row = 0; row < m_rows; row++) {
00157 for (col = 0; col < m_cols; col++) {
00158 inv.set (row, col, adjoint.get (row, col) / determ);
00159 }
00160 }
00161 } else {
00162 inv.set (0, 0, 1.0 / get (0, 0));
00163 }
00164
00165 return inv;
00166 }
00167
00168 Matrix Matrix::inverseGaussianElimination () const
00169 {
00170
00171
00172 int row, col, rowFrom, rowTo;
00173
00174
00175 QVector<int> rowIndexes (rows ());
00176 for (row = 0; row < rows (); row++) {
00177 rowIndexes [row] = row;
00178 }
00179
00180
00181
00182 Matrix working (rows (), 2 * cols ());
00183 for (row = 0; row < rows (); row++) {
00184 for (col = 0; col < cols (); col++) {
00185 double value = get (row, col);
00186 working.set (row, col, value);
00187 }
00188 }
00189
00190
00191 for (int row1 = 0; row1 < rows (); row1++) {
00192 for (int row2 = row1 + 1; row2 < rows (); row2++) {
00193 if (working.leadingZeros (row1) > working.leadingZeros (row2)) {
00194 working.switchRows (row1, row2);
00195
00196
00197 int row1Index = rowIndexes [row1];
00198 int row2Index = rowIndexes [row2];
00199 rowIndexes [row1] = row2Index;
00200 rowIndexes [row2] = row1Index;
00201 }
00202 }
00203 }
00204
00205
00206 for (row = 0; row < rows (); row++) {
00207 for (col = cols (); col < 2 * cols (); col++) {
00208 int colIdentitySubmatrix = col - cols ();
00209 double value = (row == colIdentitySubmatrix) ? 1 : 0;
00210 working.set (row, col, value);
00211 }
00212 }
00213 \
00214
00215 for (rowFrom = 0; rowFrom < rows (); rowFrom++) {
00216 int colFirstWithNonZero = rowFrom;
00217
00218
00219
00220 if (working.get (rowFrom, colFirstWithNonZero) == 0) {
00221
00222
00223 break;
00224 }
00225
00226
00227 working.normalizeRow (rowFrom, colFirstWithNonZero);
00228
00229
00230 for (rowTo = rowFrom + 1; rowTo < rows (); rowTo++) {
00231
00232 if (working.get (rowTo, colFirstWithNonZero) != 0) {
00233
00234
00235 double factor = -1.0 * working.get (rowTo, colFirstWithNonZero) / working.get (rowFrom, colFirstWithNonZero);
00236 working.addRowToAnotherWithScaling (rowFrom, rowTo, factor);
00237 }
00238 }
00239 }
00240
00241
00242 for (rowFrom = rows () - 1; rowFrom >= 0; rowFrom--) {
00243 int colFirstWithNonZero = rowFrom;
00244
00245
00246 working.normalizeRow (rowFrom, colFirstWithNonZero);
00247
00248
00249 for (rowTo = rowFrom - 1; rowTo >= 0; rowTo--) {
00250
00251 if (working.get (rowTo, colFirstWithNonZero) != 0) {
00252
00253
00254 double factor = -1.0 * working.get (rowTo, colFirstWithNonZero) / working.get (rowFrom, colFirstWithNonZero);
00255 working.addRowToAnotherWithScaling (rowFrom, rowTo, factor);
00256 }
00257 }
00258 }
00259
00260
00261 Matrix inv (working.rows ());
00262
00263 for (row = 0; row < working.rows (); row++) {
00264
00265 int rowBeforeSort = rowIndexes [row];
00266
00267 for (col = 0; col < working.rows (); col++) {
00268 int colRightHalf = col + inv.cols ();
00269 double value = working.get (rowBeforeSort, colRightHalf);
00270 inv.set (row, col, value);
00271 }
00272 }
00273
00274 return inv;
00275 }
00276
00277 unsigned int Matrix::leadingZeros (int row) const
00278 {
00279 unsigned int sum = 0;
00280 for (int col = 0; col < cols (); col++) {
00281 if (get (row, col) != 0) {
00282 ++sum;
00283 }
00284 }
00285
00286 return sum;
00287 }
00288
00289 Matrix Matrix::minorReduced (int rowOmit, int colOmit) const
00290 {
00291 ENGAUGE_ASSERT (m_rows == m_cols);
00292
00293 Matrix outMinor (m_rows - 1);
00294 int rowMinor = 0;
00295 for (int row = 0; row < m_rows; row++) {
00296
00297 if (row != rowOmit) {
00298
00299 int colMinor = 0;
00300 for (int col = 0; col < m_cols; col++) {
00301
00302 if (col != colOmit) {
00303
00304 outMinor.set (rowMinor, colMinor, get (row, col));
00305 ++colMinor;
00306 }
00307 }
00308 ++rowMinor;
00309 }
00310 }
00311
00312 return outMinor;
00313 }
00314
00315 void Matrix::normalizeRow (int rowToNormalize,
00316 int colToNormalize)
00317 {
00318 ENGAUGE_ASSERT (get (rowToNormalize, colToNormalize) != 0);
00319
00320 double factor = 1.0 / get (rowToNormalize, colToNormalize);
00321 for (int col = 0; col < cols (); col++) {
00322 double value = get (rowToNormalize, col);
00323 set (rowToNormalize, col, factor * value);
00324 }
00325 }
00326
00327 Matrix Matrix::operator* (const Matrix &other) const
00328 {
00329 ENGAUGE_ASSERT (m_cols == other.rows ());
00330
00331 Matrix out (m_rows, other.cols ());
00332
00333 for (int row = 0; row < m_rows; row++) {
00334 for (int col = 0; col < other.cols (); col++) {
00335 double sum = 0;
00336 for (int index = 0; index < m_cols; index++) {
00337 sum += get (row, index) * other.get (index, col);
00338 }
00339 out.set (row, col, sum);
00340 }
00341 }
00342
00343 return out;
00344 }
00345
00346 QVector<double> Matrix::operator* (const QVector<double> other) const
00347 {
00348 ENGAUGE_ASSERT (m_cols == other.size ());
00349
00350 QVector<double> out;
00351 out.resize (m_rows);
00352 for (int row = 0; row < m_rows; row++) {
00353 double sum = 0;
00354 for (int col = 0; col < m_cols; col++) {
00355 sum += get (row, col) * other [col];
00356 }
00357
00358 out [row] = sum;
00359 }
00360
00361 return out;
00362 }
00363
00364 int Matrix::rows () const
00365 {
00366 return m_rows;
00367 }
00368
00369 void Matrix::set (int row, int col, double value)
00370 {
00371 m_vector [fold2dIndexes (row, col)] = value;
00372 }
00373
00374 void Matrix::switchRows (int row1,
00375 int row2)
00376 {
00377
00378 for (int col = 0; col < cols (); col++) {
00379 double temp1 = get (row1, col);
00380 double temp2 = get (row2, col);
00381 set (row2, col, temp2);
00382 set (row1, col, temp1);
00383 }
00384 }
00385
00386 QString Matrix::toString () const
00387 {
00388 QString out;
00389 QTextStream str (&out);
00390
00391 str << "(";
00392 for (int row = 0; row < rows (); row++) {
00393 if (row > 0) {
00394 str << ", ";
00395 }
00396 str << "(";
00397 for (int col = 0; col < cols (); col++) {
00398 if (col > 0) {
00399 str << ", ";
00400 }
00401 str << get (row, col);
00402 }
00403 str << ")";
00404 }
00405 str << ")";
00406
00407 return out;
00408 }
00409
00410 Matrix Matrix::transpose () const
00411 {
00412 Matrix out (m_cols, m_rows);
00413
00414 for (int row = 0; row < m_rows; row++) {
00415 for (int col = 0; col < m_cols; col++) {
00416 out.set (col, row, get (row, col));
00417 }
00418 }
00419
00420 return out;
00421 }