矩阵的知识点之多足以写成一本线性代数。
在C++中,我们把矩阵封装成类。。
程序清单:
Matrix.h//未完待续
#ifndef _MATRIX_H #define _MATRIX_H #include<iostream> #include<vector> using namespace std; template <typename T> class Matrix { public://矩阵基本运算 Matrix operator*(const Matrix<T> &matrix); //运算符重载“*”重载为点乘 Matrix operator+(const Matrix<T> &matrix); //运算符重载“+”为矩阵加法 Matrix operator-(const Matrix<T> &matrix); //运算符重载“-”为矩阵减法 Matrix operator|(const Matrix<T> &matrix); //重载运算符“|”为为矩阵增广(把b矩阵并在a矩阵后面) Matrix operator*(const double input); //重载运算符“*”常数乘法 Matrix operator*(const float input); //同上 Matrix operator*(const int input); //同上 double determinant(); //求矩阵行列式的值 double det_recursion(Matrix<T> &matrix); //求矩阵行列式的值 递归 Matrix<T> getAlgebraicCofactor(int row,int column); //获取代数余子式 bool calAdjugateMatrix(Matrix<double> *matrix); //求伴随矩阵 Matrix transpose(); //获取矩阵的转置 Matrix reduceAnIdentity(int number); //产生一个秩为number的单位矩阵 bool GaussianElimination(); //高斯消元 bool rref(); //化简矩阵(reduced row echelon form)行最简式 bool rrefmovie(); //化简矩阵(reduced row echelon form)行最简式,并打印过程 bool inverse(); //矩阵的逆 bool clear(); //清空矩阵 public://矩阵增删补改 Matrix(); Matrix(int rows,int columns); //产生一个行数为rows,列数为columnns的矩阵 virtual ~Matrix(); int getRows()const; //获取行数 int getColumns()const; //获取列数 bool getSpecifiedElem(int row,int column, T *elem)const; //获取第i行第j个元素 bool getSpecifiedRow(int index,vector<T> *vec)const; //获取第index行元素 bool getSpecifiedColumn(int index,vector<T> *vec)const; //获取第index列元素 bool setSpecifiedElem(int row,int column,T value); //设置第row行第column个元素的值为value bool addOneRowToBack(vector<T> &vec); //往最底插入一行 bool addOneColumToBack(vector<T> &vec); //往最后插入一列 bool exchangeRows(int index1st,int index2nd); //交换矩阵的两行 bool relevantCheck(vector<T> &vecA,vector<T> &vecB); //两个向量的相关性检查 bool invertibleCheck(); //可逆性检查 void printfAll(); //打印所有元素 private: int m_iRows; int m_iColumns; vector<vector<T>> m_vecMatrix; }; template <typename T> Matrix<T>::Matrix() { m_iRows = 0; m_iColumns =0; m_vecMatrix.clear(); } template <typename T> Matrix<T>::Matrix(int rows,int columns) { vector<T> tempVec; m_iRows = 0; m_iColumns =0; m_vecMatrix.clear(); //reduce a vector for adding for(int i=0;i<columns;i++) { tempVec.push_back(0); } //puduce a zero matirx for(int i=0;i<rows;i++) { addOneRowToBack(tempVec); } } template <typename T> Matrix<T>::~Matrix() { m_vecMatrix.clear(); } template <typename T> int Matrix<T>::getRows()const { return m_iRows; } template <typename T> int Matrix<T>::getColumns()const { return m_iColumns; } template <typename T> bool Matrix<T>::getSpecifiedElem(int row,int column, T *elem)const { if(row>=m_iRows || column>=m_iColumns) { return false; } *elem = this->m_vecMatrix[row][column]; return true; } template <typename T> bool Matrix<T>::getSpecifiedRow(int index,vector<T> *vec)const { if(index > m_iRows) { return false; } vec->clear(); for(int i=0;i<m_iColumns;i++) { vec->push_back(m_vecMatrix[index][i]); } return true; } template <typename T> bool Matrix<T>::getSpecifiedColumn(int index,vector<T> *vec)const//获取第index列元素 { if(index > m_iColumns) { return false; } vec->clear(); for(int i=0;i<m_iRows;i++) { vec->push_back(m_vecMatrix[i][index]); } return true; } template <typename T> bool Matrix<T>::setSpecifiedElem(int row,int column,T value) //设置第row行第column个元素的值为value { if(row > m_iRows-1 ||column > m_iColumns-1) { return false; } m_vecMatrix[row][column]=value; return true; } template <typename T> bool Matrix<T>::addOneRowToBack(vector<T> &vec) { /*Matrix has had datas*/ if(m_iColumns !=0) { if(vec.size() != m_iColumns)//input data columns not equal matrix columns { cout<<"addOneRowToBack(): columns not equal"<<endl; return false ; } } /*Adding vec to m_vecDatas*/ m_vecMatrix.push_back(vec); m_iRows ++; m_iColumns = vec.size(); return true; } template <typename T> bool Matrix<T>::addOneColumToBack(vector<T> &vec) { /*Matrix has had datas*/ if(m_iRows !=0) { if(vec.size() != m_iRows)//input data columns not equal matrix columns { cout<<"addOneColumsToBack(): rows not equal"<<endl; return false ; } } else { vector<T> tempVec; m_iRows = vec.size(); for(int i=0;i<m_iRows;i++) { m_vecMatrix.push_back(tempVec); } } /*Adding vec to m_vecDatas*/ for(int i=0;i<m_iRows;i++) { m_vecMatrix[i].push_back(vec[i]); } m_iColumns ++; m_iRows = vec.size(); return true; } template <typename T> bool Matrix<T>::exchangeRows(int index1st,int index2nd) { /*input leagality check*/ if(index1st >= m_iRows || index2nd >= m_iRows) { cout<<"exchangeRows():index illeagel"; return false; } /*Exchange*/ Matrix<T> outputMatrix = *this; vector<T> tempVec = outputMatrix.m_vecMatrix[index1st]; outputMatrix.m_vecMatrix[index1st] = outputMatrix.m_vecMatrix[index2nd]; outputMatrix.m_vecMatrix[index2nd] = tempVec; *this = outputMatrix; return true; } template <typename T> Matrix<T> Matrix<T>::reduceAnIdentity(int number) { Matrix<T> outputMatrix(number,number); if(number <= 0) { return outputMatrix; } //change matrixdata[i][i] to 1 for(int i=0;i<number;i++) { outputMatrix.m_vecMatrix[i][i] = (T)1; } return outputMatrix; } template <typename T> bool Matrix<T>::GaussianElimination() { Matrix<T> outputMatrix = *this; /*Gaussian elmiation*/ for(int k=0;k<outputMatrix.m_iRows;k++) { /*if all the pivot have been found*/ if(k>=m_iColumns) { break; } /*exchange rows downward to find the row's pivot*/ for(int i=k+1;i<outputMatrix.m_iRows;i++) { /*pivot is non-zero*/ if(outputMatrix.m_vecMatrix[k][k] != 0) { //T temp = outputMatrix.m_vecMatrix[0][0]; break; } else { if(i < outputMatrix.m_iRows) { outputMatrix.exchangeRows(k,i); } } } /*if there is no pivot in this row*/ if(outputMatrix.m_vecMatrix[k][k] == 0) { break; } /*elimination:The rows below pivot row subtract times of pivot row*/ for(int i=k+1;i<outputMatrix.m_iRows;i++) { double RowsfirstData = outputMatrix.m_vecMatrix[i][k]/outputMatrix.m_vecMatrix[k][k];//Save the first data of next(k+1) rows if(RowsfirstData != 0) { outputMatrix.m_vecMatrix[i][k]=0; for(int j=k+1;j<outputMatrix.m_iColumns;j++) { outputMatrix.m_vecMatrix[i][j] -= RowsfirstData*outputMatrix.m_vecMatrix[k][j] ; } } } } *this = outputMatrix; return true; } template <typename T> bool Matrix<T>::rref() { Matrix<T> outputMatrix = *this; int rank=0;//the rank of the matrix, how many columns's pivot will it has(-1) /*Gaussian elmiation*/ for(int k=0;k<outputMatrix.m_iRows;k++) { /*if all the pivot elem have been found*/ if(k>=m_iColumns) { break; } /*exchange rows downward to find the pivot row*/ for(int i=k+1;i<outputMatrix.m_iRows;i++) { /*pivot is non-zero*/ if(outputMatrix.m_vecMatrix[k][k] != 0) { //T temp = outputMatrix.m_vecMatrix[0][0]; rank++; break; } else { if(i < outputMatrix.m_iRows) { outputMatrix.exchangeRows(k,i); } } } /*if there is no pivot in this row*/ if(outputMatrix.m_vecMatrix[k][k] == 0) { break; } /*elimination:The rows below pivot row subtract times of pivot row*/ for(int i=k+1;i<outputMatrix.m_iRows;i++) { double RowsfirstData = outputMatrix.m_vecMatrix[i][k]/outputMatrix.m_vecMatrix[k][k];//Save the first data of next(k+1) rows if(RowsfirstData != 0) { outputMatrix.m_vecMatrix[i][k]=0; for(int j=k+1;j<outputMatrix.m_iColumns;j++) { outputMatrix.m_vecMatrix[i][j] -= RowsfirstData*outputMatrix.m_vecMatrix[k][j] ; } } } } /*normalizing:set all pivots to 1*/ for(int i=0;i<outputMatrix.m_iRows;i++) { for(int j=0;j<outputMatrix.m_iColumns;j++) { if(outputMatrix.m_vecMatrix[i][j] !=0 )//pivot has been foound { double pivot = outputMatrix.m_vecMatrix[i][j];//get pivot for(int k=i;k<outputMatrix.m_iColumns;k++) { outputMatrix.m_vecMatrix[i][k] /=pivot; } break; } } } /*Back substitution*/ for(int i = rank;i>=1;i--) { /*find a first non-zero elem (It is pivot)*/ for(int j=0;j<outputMatrix.m_iColumns;j++) { double times=0; if(outputMatrix.m_vecMatrix[i][j] !=0)//pivot found { for(int l=i-1;l>=0;l--) { times = outputMatrix.m_vecMatrix[l][j]/outputMatrix.m_vecMatrix[i][j]; for(int k=j;k<outputMatrix.m_iColumns;k++)//tims of this row subtract by each columns in upon row { outputMatrix.m_vecMatrix[l][k] -= times*outputMatrix.m_vecMatrix[i][k]; } } break; } } } *this = outputMatrix; return true; } template <typename T> bool Matrix<T>::rrefmovie() { Matrix<T> outputMatrix = *this; int rank=0;//the rank of the matrix, how many columns's pivot will it has(-1) /*Gauss elmiation*/ cout<<"Gauss elimination:"<<endl; outputMatrix.printfAll(); for(int k=0;k<outputMatrix.m_iRows;k++) { /*If all the pivot elem have been found*/ if(k>=m_iColumns) { break; } /*Exchange rows downward to find the pivot row*/ for(int i=k+1;i<outputMatrix.m_iRows;i++) { /*Pivot is non-zero*/ if(outputMatrix.m_vecMatrix[k][k] != 0) { rank++; break; } else { if(i < outputMatrix.m_iRows) { outputMatrix.exchangeRows(k,i); } } if(k!=i) { cout<<"row"<<k+1<<" exchange row"<<i<<endl;//Debug outputMatrix.printfAll(); } } /*If there is no pivot in this row*/ if(outputMatrix.m_vecMatrix[k][k] == 0) { break; } /*Elimination:The rows below pivot row subtract times of pivot row*/ for(int i=k+1;i<outputMatrix.m_iRows;i++) { double RowsfirstData = outputMatrix.m_vecMatrix[i][k]/outputMatrix.m_vecMatrix[k][k];//Save the first data of next(k+1) rows if(RowsfirstData != 0) { outputMatrix.m_vecMatrix[i][k]=0; for(int j=k+1;j<outputMatrix.m_iColumns;j++) { outputMatrix.m_vecMatrix[i][j] -= RowsfirstData*outputMatrix.m_vecMatrix[k][j] ; } } cout<<"row"<<i+1<<" - "<<RowsfirstData<<"*"<<"row"<<k+1<<endl;//Debug outputMatrix.printfAll(); } } /*Normalizing:set all rows pivot to 1*/ for(int i=0;i<outputMatrix.m_iRows;i++) { for(int j=0;j<outputMatrix.m_iColumns;j++) { if(outputMatrix.m_vecMatrix[i][j] !=0 )//pivot has been foound { double pivot = outputMatrix.m_vecMatrix[i][j];//get pivot for(int k=i;k<outputMatrix.m_iColumns;k++) { outputMatrix.m_vecMatrix[i][k] /=pivot; } cout<<"row"<<i+1<<" / "<<pivot<<endl;//Debug outputMatrix.printfAll();//Debug break; } } } /*Back substitution*/ cout<<"Back substitution:"<<endl; for(int i = rank;i>=1;i--) { /*find a first non-zero elem (It is pivot)*/ for(int j=0;j<outputMatrix.m_iColumns;j++) { double times=0; if(outputMatrix.m_vecMatrix[i][j] !=0)//pivot found { for(int l=i-1;l>=0;l--) { times = outputMatrix.m_vecMatrix[l][j]/outputMatrix.m_vecMatrix[i][j]; for(int k=j;k<outputMatrix.m_iColumns;k++)//tims of this row subtract by each columns in upon row { outputMatrix.m_vecMatrix[l][k] -= times*outputMatrix.m_vecMatrix[i][k]; } cout<<"row"<<l+1<<" - "<<times<<"*"<<"row"<<i+1<<endl; outputMatrix.printfAll(); } break; } } } *this = outputMatrix; return true; } template <typename T> void Matrix<T>::printfAll() { for(int j=0;j<m_iRows;j++) { cout<<"| "; for(int i=0;i<m_iColumns;i++) { cout<<m_vecMatrix[j][i]<<" "; } cout<<"|"<<endl; } cout<<" "<<endl; } template <typename T> Matrix<T> Matrix<T>::operator+(const Matrix<T> &matrix) //运算符重载“+”为矩阵加法 { /*matrix leagality check*/ if(this->m_iRows != matrix.getRows() || this->m_iColumns != matrix.getColumns()) { return *this; } Matrix<T> outputMatrix; vector<T> tempVec; for(int i=0;i<this->m_iRows;i++) { tempVec.clear(); for(int j=0;j<this->m_iColumns;j++) { tempVec.push_back(this->m_vecMatrix[i][j] + matrix.m_vecMatrix[i][j]); } outputMatrix.addOneRowToBack(tempVec); } return outputMatrix; } template <typename T> Matrix<T> Matrix<T>::operator-(const Matrix<T> &matrix) //运算符重载“-”为矩阵减法 { /*matrix leagality check*/ if(this->m_iRows != matrix.getRows() || this->m_iColumns != matrix.getColumns()) { return *this; } Matrix<T> outputMatrix; vector<T> tempVec; for(int i=0;i<this->m_iRows;i++) { tempVec.clear(); for(int j=0;j<this->m_iColumns;j++) { tempVec.push_back(this->m_vecMatrix[i][j] - matrix.m_vecMatrix[i][j]); } outputMatrix.addOneRowToBack(tempVec); } return outputMatrix; } template <typename T> Matrix<T> Matrix<T>::operator|(const Matrix<T> &matrix) { /*input size check*/ if(this->m_iRows != matrix.getRows() ) { return *this; } /**/ Matrix<T> outputMatrix = *this; for(int i=0;i<m_iColumns;i++) { vector<T> tempVec ; if(!matrix.getSpecifiedColumn(i,&tempVec)) { return outputMatrix; } outputMatrix.addOneColumToBack(tempVec); } return outputMatrix; } template <typename T> Matrix<T> Matrix<T>::operator*(const Matrix<T> &matrix) //运算符重载*重载为点乘 { /*matrix leagality check*/ if(this->m_iColumns != matrix.getRows()) { cout<<"operator*():input ileagal"<<endl; return *this; } /*Caculate point multiply*/ //function 1st Matrix<T> outputMatrix; vector<T> tempVec; T tempData; for(int j=0;j<m_iRows;j++) { for(int k=0;k<matrix.m_iColumns;k++) { tempData =0; for(int i=0;i<m_iColumns;i++) { tempData += this->m_vecMatrix[j][i] * matrix.m_vecMatrix[i][k]; } tempVec.push_back(tempData); } outputMatrix.addOneRowToBack(tempVec); tempVec.clear(); //clear for next rows adding } //function 2nd //vector<T> tempVecRow; //vector<T> tempVecColumn; //vector<T> tempVecResult; //int temp; //for(int k=0;k<m_iRows;k++) //{ // for(int i=0;i<matrix.getColumns();i++) // { // this->getSpecifiedRow(k,&tempVecRow);//get multiplier row // matrix.getSpecifiedColumn(i,&tempVecColumn);//get multiplicand column // temp=0; // for(int j=0;j<(int)tempVecRow.size();j++) // { // temp += tempVecRow[j]*tempVecColumn[j]; // } // tempVecResult.push_back(temp); // } // outputMatrix.addOneRowToBack(tempVecResult); // tempVecResult.clear(); //clear for next rows adding //} return outputMatrix; } /*常数乘法*/ template <typename T> Matrix<T> Matrix<T>::operator*(const double input) //重载运算符“*”常数乘法 { Matrix<T> outputMatrix = *this; for(int i=0;i<m_iRows;i++) { for(int j=0;j<m_iColumns;j++) { outputMatrix.m_vecMatrix[i][j] *= input; } } return outputMatrix; } template <typename T> Matrix<T> Matrix<T>::operator*(const float input) { Matrix<T> outputMatrix = *this; for(int i=0;i<m_iRows;i++) { for(int j=0;j<m_iColumns;j++) { outputMatrix.m_vecMatrix[i][j] *= input; } } return outputMatrix; } template <typename T> Matrix<T> Matrix<T>::operator*(const int input) { Matrix<T> outputMatrix = *this; for(int i=0;i<m_iRows;i++) { for(int j=0;j<m_iColumns;j++) { outputMatrix.m_vecMatrix[i][j] *= input; } } return outputMatrix; } template <typename T> double Matrix<T>::determinant() { //not square if(m_iColumns != m_iRows) return 0; /*function 1st:recursion*/ //double result = det_recursion(*this); /*function 2nd*/ //1:Gaussian elimination //2:Diagonal multipication Matrix<T> tempMatrix = *this; tempMatrix.GaussianElimination(); double result=1; for(int i=0;i<m_iRows;i++) { result*=tempMatrix.m_vecMatrix[i][i]; } return result; } template <typename T> double Matrix<T>::det_recursion(Matrix<T> &matrix) { double result=0; if(matrix.m_iRows!=2) { for(int i=0;i<matrix.m_iColumns;i++) { Matrix<T> tempMatrix = matrix.getAlgebraicCofactor(i,0); if(i%2==0)//even { result+=matrix.m_vecMatrix[i][0]*det_recursion(tempMatrix); } else//old { result-=matrix.m_vecMatrix[i][0]*det_recursion(tempMatrix); } } return result; } else { return (double)matrix.m_vecMatrix[0][0]*matrix.m_vecMatrix[1][1] - (double)matrix.m_vecMatrix[0][1]*matrix.m_vecMatrix[1][0]; } } template <typename T> Matrix<T> Matrix<T>::getAlgebraicCofactor(int row,int column) { if(row<0 || column<0 || row >m_iColumns-1 || column>m_iColumns-1) { return *this; } Matrix<T> outputMatrix(this->m_iRows-1,this->m_iColumns-1); for(int i=0;i<this->m_iRows;i++) { for(int j=0;j<this->m_iColumns;j++) { if(i<row ) { if(j<column) outputMatrix.m_vecMatrix[i][j] = this->m_vecMatrix[i][j]; else if(j>column) outputMatrix.m_vecMatrix[i][j-1] = this->m_vecMatrix[i][j]; } else if(i>row) { if(j<column) outputMatrix.m_vecMatrix[i-1][j] = this->m_vecMatrix[i][j]; else if(j>column) outputMatrix.m_vecMatrix[i-1][j-1] = this->m_vecMatrix[i][j]; } } } return outputMatrix; } template <typename T> bool Matrix<T>::calAdjugateMatrix(Matrix<double> *matrix) { if(determinant() ==0 ) return false ; Matrix<double> outputMatrix(m_iRows,m_iColumns); for(int i=0;i<m_iRows;i++) { for(int j=0;j<m_iColumns;j++) { outputMatrix.m_vecMatrix[i][j] = getAlgebraicCofactor(i,j).determinant(); } } *matrix = outputMatrix; return true; } template <typename T> Matrix<T> Matrix<T>::transpose() { Matrix<T> tempMatrix; vector<T> tempVec; /*get transpose*/ for(int i=0;i<m_iColumns;i++) { this->getSpecifiedRow(i,&tempVec); tempMatrix.addOneColumToBack(tempVec); } /*swap rows number and columns number*/ m_vecMatrix = tempMatrix.m_vecMatrix; int temp = m_iColumns; m_iColumns = m_iRows; m_iRows = temp; return tempMatrix; } template <typename T> bool Matrix<T>::relevantCheck(vector<T> &vecA,vector<T> &vecB) { if(vecA.size()!=vecB.size()) { return false; } int length = vecA.size(); T timesofA,timesofB; int i; for(i=0;i<length;i++) { if(vecA[i]!=0 || vecB[i]!=0) { timesofA=vecA[i]; timesofB=vecB[i]; break; } if(i == length-1) { return true; } } i+=1; for(i;i<length;i++) { if(timesofA*vecB[i] != timesofB*vecA[i]) return false; } return true; } template <typename T> bool Matrix<T>::invertibleCheck() { /*rows relevant check*/ for(int i=0;i<m_iRows-1;i++) { vector<T> vecA; vector<T> vecB; for(int j=i+1;j<m_iRows;j++) { getSpecifiedRow(i,&vecA); getSpecifiedRow(j,&vecB); if(relevantCheck(vecA,vecB)) { return false; } } } /*columns relevant check*/ for(int i=0;i<m_iRows-1;i++) { vector<T> vecA; vector<T> vecB; for(int j=i+1;j<m_iRows;j++) { getSpecifiedColumn(i,&vecA); getSpecifiedColumn(j,&vecB); if(relevantCheck(vecA,vecB)) { return false; } } } return true; } template <typename T> bool Matrix<T>::inverse() { /*Determinant equart 0,there is no inverse.*/ if(this->determinant() ==0 ) { return false; } /*Make an Identity*/ Matrix<T> I =this->reduceAnIdentity( this->m_iRows ); /*Augmented by I and do rref*/ Matrix<T> augmentationMatrix=*this|I; augmentationMatrix.rref(); //augmentationMatrix.rrefmovie(); Matrix<T> outputMatrix; vector<T> tempVec; for(int i=0;i<m_iColumns;i++) { augmentationMatrix.getSpecifiedColumn(m_iRows+i,&tempVec); outputMatrix.addOneColumToBack(tempVec); } *this = outputMatrix; return true; } template <typename T> bool Matrix<T>::clear() { m_vecMatrix.clear(); return true; } #endif
分组解析:
基本成员解析:
- 矩阵 -【2】 矩阵生成:http://www.cnblogs.com/HongYi-Liang/p/7275278.html
算术运算:
- 矩阵 - 【3】矩阵加减:http://www.cnblogs.com/HongYi-Liang/p/7287403.html
- 矩阵 - 【4】矩阵点乘:http://www.cnblogs.com/HongYi-Liang/p/7287324.html
- 矩阵 - 叉乘:
其他运算:
- 矩阵 - 矩阵转置:http://www.cnblogs.com/HongYi-Liang/p/7287495.html
- 矩阵 - 【5】矩阵化简:http://www.cnblogs.com/HongYi-Liang/p/7464850.html
- 矩阵 - 行列式:
- 矩阵 - 逆矩阵: