00001
00002
00003
00004
00005 #include <iostream>
00006 #include <fstream>
00007 #include <cmath>
00008 #include "matrix.h"
00009 #include "utils.h"
00010 using namespace std;
00011
00012
00013
00014
00015 Matrix& Matrix::RightAppendIdentity()
00016 {
00017 Matrix temp;
00018 temp = Matrix(0.0, m_nRows, (2 * m_nCols));
00019
00020 for(int i=0; i<m_nRows; ++i) {
00021 for(int j=0; j<m_nCols; ++j) {
00022 temp.m_pData[i][j] = m_pData[i][j];
00023 }
00024 }
00025 for(int q=0; q<m_nRows; ++q) {
00026 temp.m_pData[q][m_nCols+q] = 1;
00027 }
00028
00029 *this = temp;
00030
00031 return *this;
00032 }
00033
00034
00035
00036
00037 Matrix& Matrix::LeftRemoveIdentity()
00038 {
00039 Matrix temp;
00040 temp = Matrix(0.0, m_nRows, (m_nCols / 2));
00041
00042 for(int i=0; i<m_nRows; ++i) {
00043 for(int j=0; j<m_nCols/2; ++j) {
00044 temp.m_pData[i][j] = m_pData[i][m_nCols/2+j];
00045 }
00046 }
00047
00048 *this = temp;
00049
00050 return *this;
00051 }
00052
00053
00054
00055 Matrix::Matrix()
00056 {
00057 m_pData = NULL;
00058 m_nCols = m_nRows = 0;
00059 }
00060
00061
00062
00063
00064 Matrix::Matrix(double InitVal, int Rows, int Cols)
00065 {
00066 m_nRows = Rows;
00067 m_nCols = Cols;
00068
00069 m_pData = new double*[m_nRows];
00070 for(int i=0; i<m_nRows; ++i) {
00071 m_pData[i] = new double[m_nCols];
00072 for(int j=0; j<m_nCols; ++j) {
00073 m_pData[i][j] = InitVal;
00074 }
00075 }
00076 }
00077
00078
00079
00080
00081 Matrix::Matrix(double* Data, int Rows, int Cols)
00082 {
00083 m_nRows = Rows;
00084 m_nCols = Cols;
00085
00086 m_pData = new double*[m_nRows];
00087 for(int i=0; i<m_nRows; ++i) {
00088 m_pData[i] = new double[m_nCols];
00089 for(int j=0; j<m_nCols; ++j) {
00090 m_pData[i][j] = Data[i*m_nRows+j];
00091 }
00092 }
00093 }
00094
00095
00096
00097
00098 Matrix::Matrix(double** Data, int Rows, int Cols)
00099 {
00100 m_nRows = Rows;
00101 m_nCols = Cols;
00102
00103 m_pData = new double*[m_nRows];
00104 for(int i=0; i<m_nRows; ++i) {
00105 m_pData[i] = new double[m_nCols];
00106 for(int j=0; j<m_nCols; ++j) {
00107 m_pData[i][j] = Data[i][j];
00108 }
00109 }
00110 }
00111
00112
00113 Matrix::Matrix(const Matrix& obj)
00114 {
00115 m_nRows = obj.m_nRows;
00116 m_nCols = obj.m_nCols;
00117
00118 m_pData = new double*[m_nRows];
00119 for(int i=0; i<m_nRows; ++i) {
00120 m_pData[i] = new double[m_nCols];
00121 for(int j=0; j<m_nCols; ++j) {
00122 m_pData[i][j] = obj.m_pData[i][j];
00123 }
00124 }
00125 }
00126
00127
00128 Matrix::~Matrix()
00129 {
00130 for(int i=0; i<m_nRows; ++i) {
00131 delete[] m_pData[i];
00132 m_pData[i] = NULL;
00133 }
00134 delete[] m_pData;
00135 m_pData = NULL;
00136 }
00137
00138
00139 Matrix& Matrix::operator +(const Matrix& obj) const
00140 {
00141 if((m_nRows != obj.m_nRows) || (m_nCols != obj.m_nCols))
00142 ErrorMsg("mismatched matrices in addition", true);
00143
00144 static Matrix temp;
00145 temp = Matrix(0.0, m_nRows, m_nCols);
00146
00147 for(int i=0; i<m_nRows; ++i) {
00148 for(int j=0; j<m_nCols; ++j) {
00149 temp.m_pData[i][j] = m_pData[i][j] + obj.m_pData[i][j];
00150 }
00151 }
00152
00153 return temp;
00154 }
00155
00156
00157 Matrix& Matrix::operator -(const Matrix& obj) const
00158 {
00159 if((m_nRows != obj.m_nRows) || (m_nCols != obj.m_nCols))
00160 ErrorMsg("mismatched matrices in addition", true);
00161
00162 static Matrix temp;
00163 temp = Matrix(0.0, m_nRows, m_nCols);
00164
00165 for(int i=0; i<m_nRows; ++i) {
00166 for(int j=0; j<m_nCols; ++j) {
00167 temp.m_pData[i][j] = m_pData[i][j] - obj.m_pData[i][j];
00168 }
00169 }
00170
00171 return temp;
00172 }
00173
00174
00175 Matrix& Matrix::operator *(const Matrix& obj) const
00176 {
00177 if(m_nCols != obj.m_nRows)
00178 ErrorMsg("mismatched matrices in multiplication", true);
00179
00180 double sum = 0;
00181 double prod = 1;
00182
00183 static Matrix temp;
00184 temp = Matrix(0.0, m_nRows, obj.m_nCols);
00185
00186 for(int i=0; i<temp.m_nRows; ++i) {
00187 for(int j=0; j<temp.m_nCols; ++j) {
00188 sum = 0;
00189 for(int q=0; q<m_nCols; ++q) {
00190 prod = m_pData[i][q] * obj.m_pData[q][j];
00191 sum += prod;
00192 }
00193 temp.m_pData[i][j] = sum;
00194 }
00195 }
00196
00197 return temp;
00198 }
00199
00200
00201 Matrix& Matrix::operator *(const double _d) const
00202 {
00203 static Matrix temp;
00204 temp = Matrix(m_pData, m_nCols, m_nRows);
00205
00206 for(int i=0; i<temp.m_nRows; ++i) {
00207 for(int j=0; j<temp.m_nCols; ++j) {
00208 temp.m_pData[i][j] *= _d;
00209 }
00210 }
00211
00212 return temp;
00213 }
00214
00215
00216 Matrix& Matrix::operator *(const int _i) const
00217 {
00218 static Matrix temp;
00219 temp = Matrix(m_pData, m_nCols, m_nRows);
00220
00221 for(int i=0; i<temp.m_nRows; ++i) {
00222 for(int j=0; j<temp.m_nCols; ++j) {
00223 temp.m_pData[i][j] *= _i;
00224 }
00225 }
00226
00227 return temp;
00228 }
00229
00230
00231 Matrix& Matrix::operator /(const Matrix& obj) const
00232 {
00233 return(*this * obj.GetInverse());
00234 }
00235
00236
00237 Matrix& Matrix::operator /(const double _d) const
00238 {
00239 if(_d == 0) ErrorMsg("divide by zero in double division", true);
00240
00241 static Matrix temp;
00242 temp = Matrix(m_pData, m_nCols, m_nRows);
00243
00244 for(int i=0; i<temp.m_nRows; ++i) {
00245 for(int j=0; j<temp.m_nCols; ++j) {
00246 temp.m_pData[i][j] /= _d;
00247 }
00248 }
00249
00250 return temp;
00251 }
00252
00253
00254 Matrix& Matrix::operator /(const int _i) const
00255 {
00256 if(_i == 0) ErrorMsg("divide by zero in integer division", true);
00257
00258 static Matrix temp;
00259 temp = Matrix(m_pData, m_nCols, m_nRows);
00260
00261 for(int i=0; i<temp.m_nRows; ++i) {
00262 for(int j=0; j<temp.m_nCols; ++j) {
00263 temp.m_pData[i][j] /= _i;
00264 }
00265 }
00266
00267 return temp;
00268 }
00269
00270
00271 Matrix& Matrix::operator +=(const Matrix& obj)
00272 {
00273 return(*this = *this + obj);
00274 }
00275
00276
00277 Matrix& Matrix::operator -=(const Matrix& obj)
00278 {
00279 return(*this = *this - obj);
00280 }
00281
00282
00283 Matrix& Matrix::operator *=(const Matrix& obj)
00284 {
00285 return(*this = *this * obj);
00286 }
00287
00288
00289 Matrix& Matrix::operator *=(const double _d)
00290 {
00291 return(*this = *this * _d);
00292 }
00293
00294
00295 Matrix& Matrix::operator *=(const int _i)
00296 {
00297 return(*this = *this * _i);
00298 }
00299
00300
00301 Matrix& Matrix::operator /=(const Matrix& obj)
00302 {
00303 return(*this = *this / obj);
00304 }
00305
00306
00307 Matrix& Matrix::operator /=(const double _d)
00308 {
00309 return(*this = *this / _d);
00310 }
00311
00312
00313 Matrix& Matrix::operator /=(const int _i)
00314 {
00315 return(*this = *this / _i);
00316 }
00317
00318
00319 Matrix& Matrix::operator =(const Matrix& obj)
00320 {
00321 if(&obj == this) return *this;
00322
00323 if(m_nRows | m_nCols) {
00324 for(int i=0; i<m_nRows; ++i) {
00325 delete[] m_pData[i];
00326 }
00327 delete[] m_pData;
00328 }
00329
00330 m_nRows = obj.m_nRows;
00331 m_nCols = obj.m_nCols;
00332
00333 m_pData = new double*[m_nRows];
00334 for(int i=0; i<m_nRows; ++i) {
00335 m_pData[i] = new double[m_nCols];
00336 for(int j=0; j<m_nCols; ++j) {
00337 m_pData[i][j] = obj.m_pData[i][j];
00338 }
00339 }
00340
00341 return *this;
00342 }
00343
00344
00345
00346 Matrix& Matrix::operator ~() const
00347 {
00348 return GetInverse();
00349 }
00350
00351
00352 bool Matrix::operator ==(const Matrix& obj) const
00353 {
00354 if(m_nRows != obj.m_nRows) return false;
00355 if(m_nCols != obj.m_nCols) return false;
00356
00357 for(int i=0; i<m_nRows; ++i) {
00358 for(int j=0; j<m_nCols; ++j) {
00359 if(m_pData[i][j] != obj.m_pData[i][j]) return false;
00360 }
00361 }
00362
00363 return true;
00364 }
00365
00366
00367 bool Matrix::operator !=(const Matrix& obj) const
00368 {
00369 if(m_nRows != obj.m_nRows) return true;
00370 if(m_nCols != obj.m_nCols) return true;
00371
00372 for(int i=0; i<m_nRows; ++i) {
00373 for(int j=0; j<m_nCols; ++j) {
00374 if(m_pData[i][j] != obj.m_pData[i][j]) return true;
00375 }
00376 }
00377
00378 return false;
00379 }
00380
00381
00382
00383 double* Matrix::operator [](const int _i) const
00384 {
00385 if((_i >= m_nRows) || (_i < 0)) ErrorMsg("out-of-bound index", true);
00386
00387 return m_pData[_i];
00388 }
00389
00390
00391
00392 double& Matrix::operator ()(const int _i, const int _j) const
00393 {
00394 if((_i >= m_nRows) || (_j >= m_nCols)) ErrorMsg("out-of-bounds index", true);
00395 if((_i < 0) || (_j < 0)) ErrorMsg("out-of-bounds index", true);
00396
00397 return m_pData[_i][_j];
00398 }
00399
00400
00401 bool Matrix::IsIdentity() const
00402 {
00403 if(m_nCols != m_nRows) return false;
00404
00405 for(int i=0; i<m_nCols; ++i) {
00406 for(int j=0; j<m_nRows; ++j) {
00407 if(i == j) {
00408 if(m_pData[i][j] != 1) return false;
00409 }
00410 else if(m_pData[i][j] != 0) return false;
00411 }
00412 }
00413
00414 return true;
00415 }
00416
00417
00418 bool Matrix::IsEmpty() const
00419 {
00420 for(int i=0; i<m_nRows; ++i) {
00421 for(int j=0; j<m_nCols; ++j) {
00422 if(m_pData[i][j] != 0) return false;
00423 }
00424 }
00425
00426 return true;
00427 }
00428
00429
00430 double Matrix::Determinant() const
00431 {
00432 if(m_nRows != m_nCols) ErrorMsg("not a square matrix for determinant", true);
00433
00434 double sum = 0;
00435
00436 if(m_nRows == 2) {
00437 return((m_pData[0][0] * m_pData[1][1]) - (m_pData[1][0] * m_pData[0][1]));
00438 }
00439
00440 for(int q=0; q<m_nCols; ++q) {
00441 Matrix* NewMinor = GetMinorNew(0, q);
00442 sum += (pow(-1, q)*(m_pData[0][q]*NewMinor->Determinant()));
00443 delete NewMinor;
00444 }
00445
00446 return sum;
00447 }
00448
00449
00450 double Matrix::SumAll() const
00451 {
00452 double sum = 0;
00453
00454 for(int i=0; i<m_nRows; ++i) {
00455 for(int j=0; j<m_nCols; ++j) {
00456 sum += m_pData[i][j];
00457 }
00458 }
00459
00460 return sum;
00461 }
00462
00463
00464 double Matrix::SumAllSquared() const
00465 {
00466 double d = SumAll();
00467
00468 return(d * d);
00469 }
00470
00471
00472 double Matrix::SumRow(const int Row) const
00473 {
00474 double sum = 0;
00475
00476 for(int j=0; j<m_nCols; ++j) {
00477 sum += m_pData[Row][j];
00478 }
00479
00480 return sum;
00481 }
00482
00483
00484 double Matrix::SumColumn(const int Col) const
00485 {
00486 double sum = 0;
00487
00488 for(int i=0; i<m_nRows; ++i) {
00489 sum += m_pData[i][Col];
00490 }
00491
00492 return sum;
00493 }
00494
00495
00496 double Matrix::SumRowSquared(const int Row) const
00497 {
00498 double d = SumRow(Row);
00499
00500 return(d * d);
00501 }
00502
00503
00504 double Matrix::SumColumnSquared(const int Col) const
00505 {
00506 double d = SumColumn(Col);
00507
00508 return(d * d);
00509 }
00510
00511
00512 double Matrix::GetMax() const
00513 {
00514 double max = m_pData[0][0];
00515
00516 for(int i=0; i<m_nRows; ++i) {
00517 for(int j=0; j<m_nCols; ++j) {
00518 if(m_pData[i][j] > max) max = m_pData[i][j];
00519 }
00520 }
00521
00522 return max;
00523 }
00524
00525
00526 double Matrix::GetMin() const
00527 {
00528 double min = m_pData[0][0];
00529
00530 for(int i=0; i<m_nRows; ++i) {
00531 for(int j=0; j<m_nCols; ++j) {
00532 if(m_pData[i][j] < min) min = m_pData[i][j];
00533 }
00534 }
00535
00536 return min;
00537 }
00538
00539
00540 double Matrix::GetRowMax(const int Row) const
00541 {
00542 double max = m_pData[Row][0];
00543
00544 for(int j=0; j<m_nCols; ++j) {
00545 if(m_pData[Row][j] > max) max = m_pData[Row][j];
00546 }
00547
00548 return max;
00549 }
00550
00551
00552 double Matrix::GetRowMin(const int Row) const
00553 {
00554 double min = m_pData[Row][0];
00555
00556 for(int j=0; j<m_nCols; ++j) {
00557 if(m_pData[Row][j] < min) min = m_pData[Row][j];
00558 }
00559
00560 return min;
00561 }
00562
00563
00564 double Matrix::GetColumnMax(const int Col) const
00565 {
00566 double max = m_pData[0][Col];
00567
00568 for(int i=0; i<m_nRows; ++i) {
00569 if(m_pData[i][Col] > max) max = m_pData[i][Col];
00570 }
00571
00572 return max;
00573 }
00574
00575
00576 double Matrix::GetColumnMin(const int Col) const
00577 {
00578 double min = m_pData[0][Col];
00579
00580 for(int i=0; i<m_nRows; ++i) {
00581 if(m_pData[i][Col] < min) min = m_pData[i][Col];
00582 }
00583
00584 return min;
00585 }
00586
00587
00588 double Matrix::GetRange() const
00589 {
00590 double min, max;
00591
00592 GetNumericRange(min, max);
00593
00594 return(max-min);
00595 }
00596
00597
00598 double Matrix::GetRowRange(const int Row) const
00599 {
00600 double min, max;
00601
00602 GetNumericRangeOfRow(min, max, Row);
00603
00604 return(max-min);
00605 }
00606
00607
00608 double Matrix::GetColumnRange(const int Col) const
00609 {
00610 double min, max;
00611
00612 GetNumericRangeOfColumn(min, max, Col);
00613
00614 return(max-min);
00615 }
00616
00617
00618
00619 double* Matrix::GetDataOneDimen() const
00620 {
00621 double* newData;
00622
00623 newData = new double[m_nRows * m_nCols];
00624
00625 for(int i=0; i<m_nRows; ++i) {
00626 for(int j=0; j<m_nCols; ++j) {
00627 newData[(i*m_nRows)+j] = m_pData[i][j];
00628 }
00629 }
00630
00631 return newData;
00632 }
00633
00634
00635
00636 double** Matrix::GetDataTwoDimen() const
00637 {
00638 double** newData;
00639
00640 newData = new double*[m_nRows];
00641 for(int i=0; i<m_nRows; ++i) {
00642 newData[i] = new double[m_nCols];
00643 for(int j=0; j<m_nCols; ++j) {
00644 newData[i][j] = m_pData[i][j];
00645 }
00646 }
00647
00648 return newData;
00649 }
00650
00651
00652 int Matrix::GetRows() const
00653 {
00654 return m_nRows;
00655 }
00656
00657
00658 int Matrix::GetColumns() const
00659 {
00660 return m_nCols;
00661 }
00662
00663
00664 Matrix& Matrix::Clear()
00665 {
00666 for(int i=0; i<m_nRows; ++i) {
00667 for(int j=0; j<m_nCols; ++j) {
00668 m_pData[i][j] = 0;
00669 }
00670 }
00671
00672 return *this;
00673 }
00674
00675
00676 Matrix& Matrix::ClearRow(const int Row)
00677 {
00678 for(int j=0; j<m_nCols; ++j) {
00679 m_pData[Row][j] = 0;
00680 }
00681
00682 return *this;
00683 }
00684
00685
00686 Matrix& Matrix::ClearColumn(const int Col)
00687 {
00688 for(int i=0; i<m_nRows; ++i) {
00689 m_pData[i][Col] = 0;
00690 }
00691
00692 return *this;
00693 }
00694
00695
00696 Matrix& Matrix::SetValue(int Row,int Col,double _d)
00697 {
00698 if(Row<m_nRows && Col < m_nCols)
00699 m_pData[Row][Col] = _d;
00700 return *this;
00701 }
00702
00703
00704
00705 Matrix& Matrix::Fill(const double _d)
00706 {
00707 for(int i=0; i<m_nRows; ++i) {
00708 for(int j=0; j<m_nCols; ++j) {
00709 m_pData[i][j] = _d;
00710 }
00711 }
00712
00713 return *this;
00714 }
00715
00716
00717 Matrix& Matrix::FillRow(const int Row, const double _d)
00718 {
00719 for(int j=0; j<m_nCols; ++j) {
00720 m_pData[Row][j] = _d;
00721 }
00722
00723 return *this;
00724 }
00725
00726
00727 Matrix& Matrix::FillColumn(const int Col, const double _d)
00728 {
00729 for(int i=0; i<m_nRows; ++i) {
00730 m_pData[i][Col] = _d;
00731 }
00732
00733 return *this;
00734 }
00735
00736
00737 Matrix& Matrix::GetInverse() const
00738 {
00739 if(m_nRows != m_nCols) {
00740 ErrorMsg("not a square matrix for inverse", true);
00741 }
00742
00743 static Matrix temp;
00744 temp = *this;
00745
00746 temp.RightAppendIdentity();
00747 temp.RREF();
00748 temp.LeftRemoveIdentity();
00749
00750 return temp;
00751 }
00752
00753
00754 Matrix& Matrix::Invert()
00755 {
00756 *this = this->GetInverse();
00757
00758 return *this;
00759 }
00760
00761
00762 Matrix& Matrix::AddRows(const int SourceRow, const int DestRow, const double factor)
00763 {
00764 for(int j=0; j<m_nCols; ++j) {
00765 m_pData[DestRow][j] += (m_pData[SourceRow][j] * factor);
00766 }
00767
00768 return *this;
00769 }
00770
00771
00772 Matrix& Matrix::MultiplyRow(const int Row, const double _d)
00773 {
00774 for(int j=0; j<m_nCols; ++j) {
00775 m_pData[Row][j] *= _d;
00776 }
00777
00778 return *this;
00779 }
00780
00781
00782 Matrix& Matrix::DivideRow(const int Row, const double _d)
00783 {
00784 if(_d == 0) ErrorMsg("divide by zero in row division", true);
00785
00786 for(int j=0; j<m_nCols; ++j) {
00787 m_pData[Row][j] /= _d;
00788 }
00789
00790 return *this;
00791 }
00792
00793
00794 Matrix& Matrix::AddColumns(const int SourceCol, const int DestCol, const double factor)
00795 {
00796 for(int i=0; i<m_nRows; ++i) {
00797 m_pData[i][DestCol] += (m_pData[i][SourceCol] * factor);
00798 }
00799
00800 return *this;
00801 }
00802
00803
00804 Matrix& Matrix::MultiplyColumn(const int Col, const double _d)
00805 {
00806 for(int i=0; i<m_nRows; ++i) {
00807 m_pData[i][Col] *= _d;
00808 }
00809
00810 return *this;
00811 }
00812
00813
00814 Matrix& Matrix::DivideColumn(const int Col, const double _d)
00815 {
00816 if(_d == 0) ErrorMsg("divide by zero in column division", true);
00817
00818 for(int i=0; i<m_nRows; ++i) {
00819 m_pData[i][Col] /= _d;
00820 }
00821
00822 return *this;
00823 }
00824
00825
00826 Matrix& Matrix::REF()
00827 {
00828 for(int i=0; i<m_nRows; ++i) {
00829 for(int j=i+1; j<m_nRows; ++j) {
00830 AddRows(i, j, -m_pData[j][i]/m_pData[i][i]);
00831 }
00832 DivideRow(i, m_pData[i][i]);
00833 }
00834
00835 return *this;
00836 }
00837
00838
00839 Matrix& Matrix::RREF()
00840 {
00841 REF();
00842
00843 for(int i=m_nRows-1; i>=0; i--) {
00844 for(int j=i-1; j>=0; j--) {
00845 AddRows(i, j, -m_pData[j][i]/m_pData[i][i]);
00846 }
00847 DivideRow(i, m_pData[i][i]);
00848 }
00849
00850 return *this;
00851 }
00852
00853
00854 Matrix& Matrix::GetREF() const
00855 {
00856 static Matrix temp;
00857
00858 temp = *this;
00859 temp.REF();
00860
00861 return temp;
00862 }
00863
00864
00865 Matrix& Matrix::GetRREF() const
00866 {
00867 static Matrix temp;
00868
00869 temp = *this;
00870 temp.RREF();
00871
00872 return temp;
00873 }
00874
00875
00876
00877 Matrix& Matrix::GetMinor(const int RowSpot, const int ColSpot) const
00878 {
00879 static Matrix temp;
00880 temp = Matrix(0.0, m_nRows-1, m_nCols-1);
00881
00882 for(int i=0, k=0; i<m_nRows; ) {
00883 if(i == RowSpot) {
00884 ++i;
00885 continue;
00886 }
00887 for(int j=0, l=0; j<m_nCols; ) {
00888 if(j == ColSpot) {
00889 ++j;
00890 continue;
00891 }
00892 temp.m_pData[k][l] = m_pData[i][j];
00893 ++l;
00894 ++j;
00895 }
00896 ++i;
00897 ++k;
00898 }
00899
00900 return temp;
00901 }
00902
00903
00904
00905 Matrix* Matrix::GetMinorNew(const int RowSpot, const int ColSpot) const
00906 {
00907 Matrix* temp;
00908 temp = new Matrix(0.0, m_nRows-1, m_nCols-1);
00909
00910 for(int i=0, k=0; i<m_nRows; ) {
00911 if(i == RowSpot) {
00912 ++i;
00913 continue;
00914 }
00915 for(int j=0, l=0; j<m_nCols; ) {
00916 if(j == ColSpot) {
00917 ++j;
00918 continue;
00919 }
00920 temp->m_pData[k][l] = m_pData[i][j];
00921 ++l;
00922 ++j;
00923 }
00924 ++i;
00925 ++k;
00926 }
00927
00928 return temp;
00929 }
00930
00931
00932
00933 Matrix& Matrix::GetSubMatrix(const int RowSpot, const int ColSpot, const int RowLen, const int ColLen) const
00934 {
00935 static Matrix temp;
00936 temp = Matrix(0.0, RowLen, ColLen);
00937
00938 for(int i=RowSpot, k=0; i<(RowLen+RowSpot); ++i, ++k) {
00939 for(int j=ColSpot, l=0; j<(ColLen+ColSpot); ++j, ++l) {
00940 temp.m_pData[k][l] = m_pData[i][j];
00941 }
00942 }
00943
00944 return temp;
00945 }
00946
00947
00948
00949 Matrix& Matrix::SetSubMatrix(const int RowSpot, const int ColSpot, const int RowLen, const int ColLen)
00950 {
00951 Matrix temp;
00952
00953 temp = GetSubMatrix(RowSpot, ColSpot, RowLen, ColLen);
00954
00955 *this = temp;
00956
00957 return *this;
00958 }
00959
00960
00961 Matrix& Matrix::SwapRows(const int Row1, const int Row2)
00962 {
00963 double* temp;
00964
00965 temp = new double[m_nCols];
00966
00967 for(int j=0; j<m_nCols; ++j) {
00968 temp[j] = m_pData[Row1][j];
00969 m_pData[Row1][j] = m_pData[Row2][j];
00970 m_pData[Row2][j] = temp[j];
00971 }
00972
00973 delete[] temp;
00974
00975 return *this;
00976 }
00977
00978
00979 Matrix& Matrix::SwapCols(const int Col1, const int Col2)
00980 {
00981 double* temp;
00982
00983 temp = new double[m_nRows];
00984
00985 for(int i=0; i<m_nRows; ++i) {
00986 temp[i] = m_pData[i][Col1];
00987 m_pData[i][Col1] = m_pData[i][Col2];
00988 m_pData[i][Col2] = temp[i];
00989 }
00990
00991 delete[] temp;
00992
00993 return *this;
00994 }
00995
00996
00997 Matrix& Matrix::GetTransposed() const
00998 {
00999 static Matrix temp;
01000 temp = Matrix(0.0, m_nCols, m_nRows);
01001
01002 for(int i=0; i<m_nRows; ++i) {
01003 for(int j=0; j<m_nCols; ++j) {
01004 temp.m_pData[j][i] = m_pData[i][j];
01005 }
01006 }
01007
01008 return temp;
01009 }
01010
01011
01012 Matrix& Matrix::Transpose()
01013 {
01014 *this = this->GetTransposed();
01015
01016 return *this;
01017 }
01018
01019
01020 Matrix& Matrix::GetNumericRange(double& Min, double& Max) const
01021 {
01022 Min = Max = m_pData[0][0];
01023
01024 for(int i=0; i<m_nRows; ++i) {
01025 for(int j=0; j<m_nCols; ++j) {
01026 if(m_pData[i][j] < Min) Min = m_pData[i][j];
01027 if(m_pData[i][j] > Max) Max = m_pData[i][j];
01028 }
01029 }
01030
01031 return const_cast<Matrix&>(*this);
01032 }
01033
01034
01035 Matrix& Matrix::GetNumericRangeOfRow(double& Min, double& Max, const int Row) const
01036 {
01037 Min = Max = m_pData[Row][0];
01038
01039 for(int j=0; j<m_nCols; ++j) {
01040 if(m_pData[Row][j] > Max) Max = m_pData[Row][j];
01041 if(m_pData[Row][j] < Min) Min = m_pData[Row][j];
01042 }
01043
01044 return const_cast<Matrix&>(*this);
01045 }
01046
01047
01048 Matrix& Matrix::GetNumericRangeOfColumn(double& Min, double& Max, const int Col) const
01049 {
01050 Min = Max = m_pData[0][Col];
01051
01052 for(int i=0; i<m_nRows; ++i) {
01053 if(m_pData[i][Col] > Max) Max = m_pData[i][Col];
01054 if(m_pData[i][Col] < Min) Min = m_pData[i][Col];
01055 }
01056
01057 return const_cast<Matrix&>(*this);
01058 }
01059
01060
01061
01062 Matrix& Matrix::CholeskyDecomposition()
01063 {
01064 int i, j, k, size = GetRows();
01065 Real sum=0.0;
01066
01067 static Matrix temp;
01068 temp = Matrix(0.0, size, size);
01069
01070 for (i=0; i<size; i++) {
01071 for (j=i; j<size; j++) {
01072 sum = m_pData[i][j];
01073 for (k=0; k<=i-1; k++) {
01074 sum -= temp[i][k]*temp[j][k];
01075 }
01076 if (i == j) {
01077
01078
01079 temp[i][i] = std::sqrt(std::max<Real>(sum, 0.0));
01080 } else {
01081
01082
01083
01084 temp[j][i] =
01085 (sum==0.0 ? 0.0 : sum/temp[i][i]);
01086 }
01087 }
01088 }
01089 return temp;
01090 }
01091
01092
01093
01094 Matrix& Matrix::CMAR(const Matrix& obj)
01095 {
01096 if(m_nCols != obj.m_nCols) ErrorMsg("mismatched matrices in row concatenation", true);
01097
01098 Matrix temp(0.0, (m_nRows + obj.m_nRows), m_nCols);
01099
01100 for(int i=0; i<m_nRows; ++i) {
01101 for(int j=0; j<m_nCols; ++j) {
01102 temp.m_pData[i][j] = m_pData[i][j];
01103 }
01104 }
01105 for(int k=0; i<temp.m_nRows; ++i, ++k) {
01106 for(int j=0; j<m_nCols; ++j) {
01107 temp.m_pData[i][j] = obj.m_pData[k][j];
01108 }
01109 }
01110
01111 *this = temp;
01112
01113 return *this;
01114 }
01115
01116
01117
01118 Matrix& Matrix::CMAC(const Matrix& obj)
01119 {
01120 if(m_nRows != obj.m_nRows) ErrorMsg("mismatched matrices in column concatenation", true);
01121
01122 Matrix temp(0.0, m_nRows, (m_nCols + obj.m_nCols));
01123
01124 for(int i=0; i<m_nRows; ++i) {
01125 for(int j=0; j<m_nCols; ++j) {
01126 temp.m_pData[i][j] = m_pData[i][j];
01127 }
01128 for(int l=0; l<obj.m_nCols; ++l, ++j) {
01129 temp.m_pData[i][j] = obj.m_pData[i][l];
01130 }
01131 }
01132
01133 *this = temp;
01134
01135 return *this;
01136 }
01137
01138
01139
01140 Matrix& Matrix::GetCMAR(const Matrix& obj) const
01141 {
01142 static Matrix temp;
01143
01144 temp = *this;
01145 temp.CMAR(obj);
01146
01147 return temp;
01148 }
01149
01150
01151
01152 Matrix& Matrix::GetCMAC(const Matrix& obj) const
01153 {
01154 static Matrix temp;
01155
01156 temp = *this;
01157 temp.CMAC(obj);
01158
01159 return temp;
01160 }
01161
01162
01163 Matrix& Matrix::ConcatenateRow(const double* RowData)
01164 {
01165 Matrix temp(0.0, m_nRows+1, m_nCols);
01166
01167 for(int i=0; i<m_nRows; ++i) {
01168 for(int j=0; j<m_nCols; ++j) {
01169 temp.m_pData[i][j] = m_pData[i][j];
01170 }
01171 }
01172 int j;
01173 for(i=m_nRows, j=0; j<m_nCols; ++j) {
01174 temp.m_pData[i][j] = RowData[j];
01175 }
01176
01177 *this = temp;
01178
01179 return *this;
01180 }
01181
01182
01183 Matrix& Matrix::ConcatenateColumn(const double* ColumnData)
01184 {
01185
01186 Matrix temp(0.0, m_nRows, m_nCols+1);
01187
01188 for(int i=0; i<m_nRows; ++i) {
01189 for(int j=0; j<m_nCols; ++j) {
01190 temp.m_pData[i][j] = m_pData[i][j];
01191 }
01192 }
01193 int j;
01194 for(j=m_nCols, i=0; i<m_nRows; ++i) {
01195 temp.m_pData[i][j] = ColumnData[i];
01196 }
01197
01198 *this = temp;
01199
01200 return *this;
01201 }
01202
01203
01204 Matrix& Matrix::SpliceInRow(const double* RowData, const int RowSpot)
01205 {
01206 Matrix temp(0.0, m_nRows+1, m_nCols);
01207
01208 for(int i=0, k=0; i<m_nRows; ++i, ++k) {
01209 if(i == RowSpot) {
01210 for(int j=0; j<m_nCols; ++j) {
01211 temp.m_pData[i][j] = RowData[j];
01212 }
01213 ++k;
01214 }
01215 for(int j=0; j<m_nCols; ++j) {
01216 temp.m_pData[k][j] = m_pData[i][j];
01217 }
01218 }
01219
01220 *this = temp;
01221
01222 return *this;
01223 }
01224
01225
01226 Matrix& Matrix::SpliceInColumn(const double* ColumnData, const int ColumnSpot)
01227 {
01228 Matrix temp(0.0, m_nRows, m_nCols+1);
01229
01230 for(int i=0; i<m_nRows; ++i) {
01231 for(int j=0, l=0; j<m_nCols; ++j, ++l) {
01232 if(j == ColumnSpot) {
01233 temp.m_pData[i][l] = ColumnData[i];
01234 ++l;
01235 }
01236 temp.m_pData[i][l] = m_pData[i][j];
01237 }
01238 }
01239
01240 *this = temp;
01241
01242 return *this;
01243 }
01244
01245
01246 Matrix& Matrix::RemoveRow(const int Row)
01247 {
01248 Matrix temp(0.0, m_nRows-1, m_nCols);
01249
01250 for(int i=0, k=0; i<m_nRows; ++i, ++k) {
01251 if(i == Row) ++i;
01252 for(int j=0; j<m_nCols; ++j) {
01253 temp.m_pData[k][j] = m_pData[i][j];
01254 }
01255 }
01256
01257 *this = temp;
01258
01259 return *this;
01260 }
01261
01262
01263 Matrix& Matrix::RemoveColumn(const int Column)
01264 {
01265 Matrix temp(0.0, m_nRows, m_nCols-1);
01266
01267 for(int i=0; i<m_nRows; ++i) {
01268 for(int j=0, l=0; j<m_nCols; ++j, ++l) {
01269 if(j == Column) ++j;
01270 temp.m_pData[i][l] = m_pData[i][j];
01271 }
01272 }
01273
01274 *this = temp;
01275
01276 return *this;
01277 }
01278
01279
01280
01281 Matrix& Matrix::SortAscend()
01282 {
01283 double* Data = GetDataOneDimen();
01284 int length = (m_nRows * m_nCols);
01285 double temp;
01286
01287 for(int i=0; i<length; ++i) {
01288 for(int j=0; j<length-1; ++j) {
01289 if(Data[j] > Data[j+1]) {
01290 temp = Data[j];
01291 Data[j] = Data[j+1];
01292 Data[j+1] = temp;
01293 }
01294 }
01295 }
01296
01297 *this = Matrix(Data, m_nRows, m_nCols);
01298
01299 return *this;
01300 }
01301
01302
01303 Matrix& Matrix::SortDescend()
01304 {
01305 double* Data = GetDataOneDimen();
01306 int length = (m_nRows * m_nCols);
01307 double temp;
01308
01309 for(int i=0; i<length; ++i) {
01310 for(int j=0; j<length-1; ++j) {
01311 if(Data[j] < Data[j+1]) {
01312 temp = Data[j];
01313 Data[j] = Data[j+1];
01314 Data[j+1] = temp;
01315 }
01316 }
01317 }
01318
01319 *this = Matrix(Data, m_nRows, m_nCols);
01320
01321 return *this;
01322 }
01323
01324
01325 Matrix& Matrix::GetNormalized(const double Min, const double Max) const
01326 {
01327 static Matrix temp;
01328
01329 temp = *this;
01330 temp.Normalize(Min, Max);
01331
01332 return temp;
01333 }
01334
01335
01336 Matrix& Matrix::Normalize(const double Min, const double Max)
01337 {
01338 double MatMin, MatMax;
01339 double Range, R_Range;
01340
01341 GetNumericRange(MatMin, MatMax);
01342
01343 Range = MatMax - MatMin;
01344 R_Range = Max - Min;
01345
01346 for(int i=0; i<m_nRows; ++i) {
01347 for(int j=0; j<m_nCols; ++j) {
01348 m_pData[i][j] -= MatMin;
01349 m_pData[i][j] /= Range;
01350 m_pData[i][j] *= R_Range;
01351 m_pData[i][j] += Min;
01352 }
01353 }
01354
01355 return *this;
01356 }
01357
01358
01359 Matrix& Matrix::GetCovariant() const
01360 {
01361 Matrix temp;
01362
01363 temp = *this;
01364 temp.Transpose();
01365
01366 return(*this * temp);
01367 }
01368
01369
01370 Matrix& Matrix::MakeCovariant()
01371 {
01372 *this = this->GetCovariant();
01373
01374 return *this;
01375 }
01376
01377
01378
01379 void Matrix::Display() const
01380 {
01381 cout << "[";
01382 for(int i=0; i<m_nRows; ++i) {
01383 for(int j=0; j<m_nCols; ++j) {
01384 cout << "[" << m_pData[i][j] << "]";
01385 }
01386 if(m_nRows-i-1) cout << "\n ";
01387 }
01388 cout << "]\n";
01389 }
01390
01391
01392
01393
01394 void Matrix::Output(ostream& ostr) const
01395 {
01396 ostr << "[";
01397 for(int i=0; i<m_nRows; ++i) {
01398 for(int j=0; j<m_nCols; ++j) {
01399 ostr << "[" << m_pData[i][j] << "]";
01400 }
01401 if(m_nRows-i-1) cout << "\n ";
01402 }
01403 ostr << "]\n";
01404 }
01405
01406 void Matrix::Input(istream& istr)
01407 {
01408 cout << "Rows: ";
01409 m_nRows = getint(istr);
01410
01411 cout << "Columns: ";
01412 m_nCols = getint(istr);
01413
01414 cout << "Data:\n";
01415 for(int i=0; i<m_nRows; ++i) {
01416 for(int j=0; j<m_nCols; ++j) {
01417 cout << "(" << i << "," << j << "): ";
01418 m_pData[i][j] = getint(istr);
01419 }
01420 }
01421 }
01422
01423
01424 void Matrix::Read(ifstream& istr)
01425 {
01426 char str[6];
01427
01428 istr >> str >> m_nRows;
01429 istr >> str >> m_nCols;
01430 istr >> str;
01431 for(int i=0; i<m_nRows; ++i) {
01432 for(int j=0; j<m_nCols; ++j) {
01433 istr >> m_pData[i][j];
01434 }
01435 }
01436 }
01437
01438
01439 void Matrix::Write(ofstream& ostr) const
01440 {
01441 ostr << "Rows: " << m_nRows << '\n';
01442 ostr << "Cols: " << m_nCols << '\n';
01443 ostr << "Data:\n";
01444 for(int i=0; i<m_nRows; ++i) {
01445 for(int j=0; j<m_nCols; ++j) {
01446 ostr << m_pData[i][j] << '\n';
01447 }
01448 }
01449 }
01450
01451
01452 Matrix& Matrix::IdentityMatrix(int Diagonal)
01453 {
01454 static Matrix temp;
01455 temp = Matrix(0.0, Diagonal, Diagonal);
01456
01457 for(int q=0; q<Diagonal; ++q) {
01458 temp.m_pData[q][q] = 1;
01459 }
01460
01461 return temp;
01462 }
01463
01464
01465 Matrix& IdentityMatrix(int Diagonal)
01466 {
01467 return Matrix::IdentityMatrix(Diagonal);
01468 }
01469
01470
01471
01472 ostream& operator <<(ostream& ostr, const Matrix& obj)
01473 {
01474 ostr << "[";
01475 for(int i=0; i<obj.GetRows(); ++i) {
01476 for(int j=0; j<obj.GetColumns(); ++j) {
01477 ostr << "[" << obj[i][j] << "]";
01478 }
01479 if(obj.GetRows()-i-1) cout << "\n ";
01480 }
01481 ostr << "]\n";
01482
01483 return ostr;
01484 }
01485
01486