矩阵的知识点之多足以写成一本线性代数。
在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 = ;
m_iColumns =;
m_vecMatrix.clear();
} template <typename T>
Matrix<T>::Matrix(int rows,int columns)
{
vector<T> tempVec; m_iRows = ;
m_iColumns =;
m_vecMatrix.clear(); //reduce a vector for adding
for(int i=;i<columns;i++)
{
tempVec.push_back();
} //puduce a zero matirx
for(int i=;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=;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=;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- ||column > m_iColumns-)
{
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 !=)
{
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 !=)
{
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=;i<m_iRows;i++)
{
m_vecMatrix.push_back(tempVec);
}
} /*Adding vec to m_vecDatas*/
for(int i=;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 <= )
{
return outputMatrix;
} //change matrixdata[i][i] to 1
for(int i=;i<number;i++)
{
outputMatrix.m_vecMatrix[i][i] = (T);
} return outputMatrix;
} template <typename T>
bool Matrix<T>::GaussianElimination()
{
Matrix<T> outputMatrix = *this; /*Gaussian elmiation*/
for(int k=;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+;i<outputMatrix.m_iRows;i++)
{
/*pivot is non-zero*/
if(outputMatrix.m_vecMatrix[k][k] != )
{
//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] == )
{
break;
} /*elimination:The rows below pivot row subtract times of pivot row*/
for(int i=k+;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 != )
{
outputMatrix.m_vecMatrix[i][k]=;
for(int j=k+;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=;//the rank of the matrix, how many columns's pivot will it has(-1) /*Gaussian elmiation*/
for(int k=;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+;i<outputMatrix.m_iRows;i++)
{
/*pivot is non-zero*/
if(outputMatrix.m_vecMatrix[k][k] != )
{
//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] == )
{
break;
} /*elimination:The rows below pivot row subtract times of pivot row*/
for(int i=k+;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 != )
{
outputMatrix.m_vecMatrix[i][k]=;
for(int j=k+;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=;i<outputMatrix.m_iRows;i++)
{
for(int j=;j<outputMatrix.m_iColumns;j++)
{
if(outputMatrix.m_vecMatrix[i][j] != )//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>=;i--)
{
/*find a first non-zero elem (It is pivot)*/
for(int j=;j<outputMatrix.m_iColumns;j++)
{
double times=;
if(outputMatrix.m_vecMatrix[i][j] !=)//pivot found
{
for(int l=i-;l>=;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=;//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=;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+;i<outputMatrix.m_iRows;i++)
{
/*Pivot is non-zero*/
if(outputMatrix.m_vecMatrix[k][k] != )
{
rank++;
break;
}
else
{
if(i < outputMatrix.m_iRows)
{
outputMatrix.exchangeRows(k,i);
}
}
if(k!=i)
{
cout<<"row"<<k+<<" exchange row"<<i<<endl;//Debug
outputMatrix.printfAll();
}
} /*If there is no pivot in this row*/
if(outputMatrix.m_vecMatrix[k][k] == )
{
break;
} /*Elimination:The rows below pivot row subtract times of pivot row*/
for(int i=k+;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 != )
{
outputMatrix.m_vecMatrix[i][k]=;
for(int j=k+;j<outputMatrix.m_iColumns;j++)
{
outputMatrix.m_vecMatrix[i][j] -= RowsfirstData*outputMatrix.m_vecMatrix[k][j] ;
}
}
cout<<"row"<<i+<<" - "<<RowsfirstData<<"*"<<"row"<<k+<<endl;//Debug
outputMatrix.printfAll();
}
} /*Normalizing:set all rows pivot to 1*/
for(int i=;i<outputMatrix.m_iRows;i++)
{
for(int j=;j<outputMatrix.m_iColumns;j++)
{
if(outputMatrix.m_vecMatrix[i][j] != )//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+<<" / "<<pivot<<endl;//Debug
outputMatrix.printfAll();//Debug
break;
}
}
} /*Back substitution*/
cout<<"Back substitution:"<<endl;
for(int i = rank;i>=;i--)
{
/*find a first non-zero elem (It is pivot)*/
for(int j=;j<outputMatrix.m_iColumns;j++)
{
double times=;
if(outputMatrix.m_vecMatrix[i][j] !=)//pivot found
{
for(int l=i-;l>=;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+<<" - "<<times<<"*"<<"row"<<i+<<endl;
outputMatrix.printfAll();
}
break;
}
}
} *this = outputMatrix;
return true;
} template <typename T>
void Matrix<T>::printfAll()
{
for(int j=;j<m_iRows;j++)
{
cout<<"| ";
for(int i=;i<m_iColumns;i++)
{
cout<<m_vecMatrix[j][i]<<" ";
}
cout<<"|"<<endl;
}
cout<<" \r\n"<<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=;i<this->m_iRows;i++)
{
tempVec.clear();
for(int j=;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=;i<this->m_iRows;i++)
{
tempVec.clear();
for(int j=;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=;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=;j<m_iRows;j++)
{
for(int k=;k<matrix.m_iColumns;k++)
{
tempData =;
for(int i=;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=;i<m_iRows;i++)
{
for(int j=;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=;i<m_iRows;i++)
{
for(int j=;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=;i<m_iRows;i++)
{
for(int j=;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 ; /*function 1st:recursion*/
//double result = det_recursion(*this); /*function 2nd*/
//1:Gaussian elimination
//2:Diagonal multipication
Matrix<T> tempMatrix = *this;
tempMatrix.GaussianElimination(); double result=;
for(int i=;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=;
if(matrix.m_iRows!=)
{
for(int i=;i<matrix.m_iColumns;i++)
{
Matrix<T> tempMatrix = matrix.getAlgebraicCofactor(i,);
if(i%==)//even
{
result+=matrix.m_vecMatrix[i][]*det_recursion(tempMatrix);
}
else//old
{
result-=matrix.m_vecMatrix[i][]*det_recursion(tempMatrix);
}
}
return result;
}
else
{
return (double)matrix.m_vecMatrix[][]*matrix.m_vecMatrix[][] - (double)matrix.m_vecMatrix[][]*matrix.m_vecMatrix[][];
}
} template <typename T>
Matrix<T> Matrix<T>::getAlgebraicCofactor(int row,int column)
{
if(row< || column< || row >m_iColumns- || column>m_iColumns-)
{
return *this;
} Matrix<T> outputMatrix(this->m_iRows-,this->m_iColumns-); for(int i=;i<this->m_iRows;i++)
{
for(int j=;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-] = this->m_vecMatrix[i][j];
}
else 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-] = this->m_vecMatrix[i][j];
} }
}
return outputMatrix;
} template <typename T>
bool Matrix<T>::calAdjugateMatrix(Matrix<double> *matrix)
{
if(determinant() == )
return false ; Matrix<double> outputMatrix(m_iRows,m_iColumns);
for(int i=;i<m_iRows;i++)
{
for(int j=;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=;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=;i<length;i++)
{
if(vecA[i]!= || vecB[i]!=)
{
timesofA=vecA[i];
timesofB=vecB[i];
break;
}
if(i == length-)
{
return true;
}
}
i+=;
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=;i<m_iRows-;i++)
{
vector<T> vecA;
vector<T> vecB;
for(int j=i+;j<m_iRows;j++)
{
getSpecifiedRow(i,&vecA);
getSpecifiedRow(j,&vecB);
if(relevantCheck(vecA,vecB))
{
return false;
}
}
} /*columns relevant check*/
for(int i=;i<m_iRows-;i++)
{
vector<T> vecA;
vector<T> vecB;
for(int j=i+;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() == )
{
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=;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
- 矩阵 - 行列式：
- 矩阵 - 逆矩阵：