Main Page | Namespace List | Class Hierarchy | Class List | File List | Class Members | File Members

matrix.cpp

Go to the documentation of this file.
00001 //------------------------------------
00002 //   Implementation of Matrix Class
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 //appends an identity matrix to the right of the calling object
00013 //used in finding the inverse
00014 //the user probably doesn't need access to this
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 //removes an identity matrix from the left of the calling object
00035 //used in finding the inverse
00036 //the user probably doesn't need access to this
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 //standard constructor
00054 //Data = Null; rows = columns = 0
00055 Matrix::Matrix()
00056 {
00057         m_pData = NULL;
00058         m_nCols = m_nRows = 0;
00059 }
00060 
00061 //constructor with explicit values
00062 //of the matrix, rows = 'Rows', columns = 'Cols'
00063 //all spots of matrix = InitVal
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 //constructor with explicit values
00079 //of the matrix, rows = 'Rows', columns = 'Cols'
00080 //all spots of matrix are the spots of 'Data', going row-by-row
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 //constructor with explicit values
00096 //of the matrix, rows = 'Rows', columns = 'Cols'
00097 //spots in matrix correspond to the spots in the double Pointer-Pointer
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 //copy constructor
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 //destructor
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 //add two matrices
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 //subtract two matrices
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 //multiply two matrices
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 //multiply a matrix by a double scalar
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 //multiply a matrix by an int scalar
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 //divide a matrix by another matrix
00231 Matrix& Matrix::operator /(const Matrix& obj) const
00232 {
00233         return(*this * obj.GetInverse());
00234 }
00235 
00236 //divide a matrix by a double scalar
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 //divide a matrix by an int scalar
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 //add two matrices and assign
00271 Matrix& Matrix::operator +=(const Matrix& obj)
00272 {
00273         return(*this = *this + obj);
00274 }
00275 
00276 //subtract two matrices and assign
00277 Matrix& Matrix::operator -=(const Matrix& obj)
00278 {
00279         return(*this = *this - obj);
00280 }
00281 
00282 //multiply two matrices and assign
00283 Matrix& Matrix::operator *=(const Matrix& obj)
00284 {
00285         return(*this = *this * obj);
00286 }
00287 
00288 //multiply a matrix by a double scalar and assign
00289 Matrix& Matrix::operator *=(const double _d)
00290 {
00291         return(*this = *this * _d);
00292 }
00293 
00294 //multiply a matrix by an int scalar and assign
00295 Matrix& Matrix::operator *=(const int _i)
00296 {
00297         return(*this = *this * _i);
00298 }
00299 
00300 //divide a matrix by a matrix and assign
00301 Matrix& Matrix::operator /=(const Matrix& obj)
00302 {
00303         return(*this = *this / obj);
00304 }
00305 
00306 //divide a matrix by a double scalar and assign
00307 Matrix& Matrix::operator /=(const double _d)
00308 {
00309         return(*this = *this / _d);
00310 }
00311 
00312 //divide a matrix by an int scalar and assign
00313 Matrix& Matrix::operator /=(const int _i)
00314 {
00315         return(*this = *this / _i);
00316 }
00317 
00318 //assignment operator
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 //inversion operator
00345 //returns the inverse of the calling matrix
00346 Matrix& Matrix::operator ~() const
00347 {
00348         return GetInverse();
00349 }
00350 
00351 //equality operator
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 //inequality operator
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 //indexing operator
00382 //remember, a zero (0) element DOES exist in this matrix, although not in a double matrix
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 //another indexing operator, but takes both rows and columns
00391 //remember, a zero (0) element DOES exist in this matrix, although not in a double matrix
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 //returns true if the calling matrix is an identity matrix
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 //returns true if every element in the calling matrix is zero (0)
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 //finds determinant of the matrix
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 //sums all the values in the matrix
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 //sums all the values in the matrix and then squares that value
00464 double Matrix::SumAllSquared() const
00465 {
00466         double d = SumAll();
00467 
00468         return(d * d);
00469 }
00470 
00471 //sums all the values in the row 'Row' of the matrix
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 //sums all the values in the column 'Col' of the matrix
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 //sums all the values in the row 'Row' of the matrix and then squares that value
00496 double Matrix::SumRowSquared(const int Row) const
00497 {
00498         double d = SumRow(Row);
00499 
00500         return(d * d);
00501 }
00502 
00503 //sums all the values in the column 'Col' of the matrix and th squares that value
00504 double Matrix::SumColumnSquared(const int Col) const
00505 {
00506         double d = SumColumn(Col);
00507 
00508         return(d * d);
00509 }
00510 
00511 //returns the largest value in the matrix
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 //returns the smallest value in the matrix
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 //returns the largest value in row 'Row' in the matrix
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 //returns the smallest value in row 'Row' in the matrix
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 //returns the largest value in column 'Col' in the matrix
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 //returns the smallest value in column 'Col' in the matrix
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 //returns the difference of the largest and smallest values in the matrix
00588 double Matrix::GetRange() const
00589 {
00590         double min, max;
00591 
00592         GetNumericRange(min, max);
00593 
00594         return(max-min);
00595 }
00596 
00597 //returns the difference of the largest and smallest values in row 'Row' in the matrix
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 //returns the difference of the largest and smallest values in column 'Col' in the matrix
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 //returns the data in a one-dimensional array
00618 //the array is dynamically allocated
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 //returns the data in a two-dimensional array
00635 //the array is dynamically allocated
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 //returns number of rows of the matrix
00652 int Matrix::GetRows() const
00653 {
00654         return m_nRows;
00655 }
00656 
00657 //returns number of columns of the matrix
00658 int Matrix::GetColumns() const
00659 {
00660         return m_nCols;
00661 }
00662 
00663 //clears every entry in the calling matrix
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 //clears every entry in the row 'Row' from the calling matrix
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 //clears every entry in the column 'Col' from the calling matrix
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 //fills the entry in the calling matrix to '_d'
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 //fills every entry in the calling matrix to '_d'
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 //fills every entry in the row 'Row' of the calling matrix to '_d'
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 //fills every entry in the column 'Col' of the calling matrix to '_d'
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 //returns the inverse of the calling matrix
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 //turns the calling matrix into its inverse
00754 Matrix& Matrix::Invert()
00755 {
00756         *this = this->GetInverse();
00757 
00758         return *this;
00759 }
00760 
00761 //add 'SourceRow' by scale of 'factor' (default is one) to 'DestRow'
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 //multiply every element in row 'Row' by '_d'
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 //divide every element in row 'Row' by '_d'
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 //add 'SourceCol' by scale of 'factor' (default is one) to 'DestCol'
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 //multiply every element in row 'Col' by '_d'
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 //divide every element in row 'Col' by '_d'
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 //puts the matrix into row-echelon form
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 //puts the matrix into reduced row-echelon form
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 //returns a matrix that is the row-echelon form of the calling matrix
00854 Matrix& Matrix::GetREF() const
00855 {
00856         static Matrix temp;
00857 
00858         temp = *this;
00859         temp.REF();
00860 
00861         return temp;
00862 }
00863 
00864 //returns a matrix that is the reduced row-echelon form of the calling matrix
00865 Matrix& Matrix::GetRREF() const
00866 {
00867         static Matrix temp;
00868 
00869         temp = *this;
00870         temp.RREF();
00871 
00872         return temp;
00873 }
00874 
00875 //returns minor around spot Row,Col
00876 //returns a static Matrix object
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 //returns minor around spot Row,Col
00904 //returns a Pointer to a dynamically allocated Matrix object
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 //returns the submatrix starting at spot (RowSpot, ColSpot) and with lengths
00932 //of 'RowLen' rows and 'ColLen' columns
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 //changes the calling matrix into a submatrix starting at spot (RowSpot, ColSpot)
00948 //and with lengths of 'RowLen' rows and 'ColLen' columns
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 //swaps two rows, Row1 and Row2, from the calling matrix
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 //swaps two columns, Col1 and Col2, from the calling matrix
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 //returns the transposition of the calling matrix
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 //transposes the calling matrix
01012 Matrix& Matrix::Transpose()
01013 {
01014         *this = this->GetTransposed();
01015 
01016         return *this;
01017 }
01018 
01019 //sets 'Min' to the minimum and 'Max to the maximum values in the matrix
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 //sets 'Min' to the minimum and 'Max to the maximum values in row 'Row' in the matrix
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 //sets 'Min' to the minimum and 'Max to the maximum values in column 'Col' in the matrix
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 // Added from various sources -- here no test for symetry
01061 // returns an upper triangular matrix such that U*U.t() = Orignial matrix
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                 // To handle positive semi-definite matrices take the
01078                 // square root of sum if positive, else zero.
01079                 temp[i][i] = std::sqrt(std::max<Real>(sum, 0.0));
01080             } else {
01081                 // With positive semi-definite matrices is possible
01082                 // to have temp[i][i]==0.0
01083                 // In this case sum happens to be zero as well
01084                 temp[j][i] =
01085                     (sum==0.0 ? 0.0 : sum/temp[i][i]);
01086             }
01087         }
01088     }
01089         return temp;
01090 }
01091 
01092 //CMAR = Concatenate Matrix As Rows
01093 //concatenates matrix 'obj' on to the right of the calling matrix
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 //CMAC = Concatenate Matrix As Columns
01117 //concatenates matrix 'obj' on to the bottom of the calling matrix
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 //CMAR = Concatenate Matrix As Rows
01139 //returns a new matrix that is the calling object + 'obj' on the right
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 //CMAC = Concatenate Matrix As Columns
01151 //returns a new matrix that is the valling object + 'obj' on the bottom
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 //adds a row onto the right of the calling matrix
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 //adds a column onto the bottom of the calling matrix
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 //adds a row into the matrix in the spot 'RowSpot'
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 //adds a column into the matrix in the spot 'ColumnSpot'
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 //removes the specified row from the calling matrix
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 //removes the specified column from the calling matrix
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 //sorts the matrix in ascending order
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 //sorts the matrix in descending order
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 //returns a matrix that is scaled between Min and Max of the calling matrix
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 //scales the values between Min and Max
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 //returns the covariant of the calling matrix (transposed(obj) * obj)
01359 Matrix& Matrix::GetCovariant() const
01360 {
01361         Matrix temp;
01362 
01363         temp = *this;
01364         temp.Transpose();
01365 
01366         return(*this * temp);
01367 }
01368 
01369 //turns this matrix into its covariant(obj = transposed(obj) * obj)
01370 Matrix& Matrix::MakeCovariant()
01371 {
01372         *this = this->GetCovariant();
01373 
01374         return *this;
01375 }
01376 
01377 //a way to display the matrix
01378 //formatted
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 //another way to display the matrix
01392 //can output to a stream other than cout, which is the default
01393 //formatted
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 //input from a file
01424 void Matrix::Read(ifstream& istr)
01425 {
01426         char str[6];
01427 
01428         istr >> str >> m_nRows;                     //"Rows: "
01429         istr >> str >> m_nCols;                     //"Cols: "
01430         istr >> str;                                //"Data\n"
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 //output to a file
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 //returns an identity matrix of size Diagonal
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 //uses Matrix::IdentityMatrix to get an identity matrix of size Diagonal
01465 Matrix& IdentityMatrix(int Diagonal)
01466 {
01467         return Matrix::IdentityMatrix(Diagonal);
01468 }
01469 
01470 //easily display the matrix, usually to the console
01471 //formatted output
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 

Note: Generated nightly - reload for latest version
Generated on Thu Dec 22 23:12:36 2005 for terreneuve by doxygen 1.3.6