c++课程设计中做了一个矩阵类，完成了一些矩阵之间的常见运算。在求特征值和特征向量部分的代码中存在一些问题，所以求特征值和特征向量的相关代码我就不放上去了，下面是矩阵其他运算的相关代码：
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<iostream>
#include<queue>
#include<windows.h>
#include<iomanip>
using namespace std;
class matrix{
int row,col;//记录行数和列数
double *p;
public:
matrix(int r=1,int c=1);//缺省的构造函数
matrix(int r,int c,double *a);//不缺参数的构造函数
matrix(const matrix &a);//拷贝构造函数
~matrix();//析构函数
friend istream& operator>>(istream& in,matrix &x);//插入运算符重载
friend ostream& operator<<(ostream& out,matrix &x);//提取运算符重载
matrix &operator=(const matrix &a);//=运算符重载函数
matrix &operator+=(const matrix &a);//+=运算符重载函数
matrix &operator-=(const matrix &a);//-=运算符重载函数
matrix &operator*=(const matrix &a);//*=运算符重载函数
matrix operator^(int n);//计算矩阵的n次幂
double &operator()(int r,int c);//()运算符重载函数
matrix Transpose();//转置函数
matrix Transformation();//返回一个非零行首非零元为1的行简化阶梯型矩阵，不改变原矩阵
double Det();//求方阵的行列式的值
friend matrix operator+(const matrix &a,const matrix &b);//+运算符重载函数
friend matrix operator-(const matrix &a,const matrix &b);//-运算符重载函数
friend matrix operator*(const matrix &a,const matrix &b);//*运算符重载函数
int rank();//求矩阵的秩
void cal(matrix b);//计算AX=b的解
matrix unit(int n);//返回n阶单位矩阵
matrix aug(const matrix &b);//返回当前矩阵与矩阵b的增广矩阵
matrix inverse(bool &flag);//求当前矩阵的逆，如若矩阵存在逆元则flag为1，反之为0
//以下三行为求矩阵特征值和特征向量的函数
friend void QR(matrix &A,matrix &Q,matrix &R);
friend void Get_Eig(const matrix &A,matrix &Vec,matrix &EigVal);
void Eig();
};
matrix::matrix(int r,int c)
{
row=r;col=c;//存储矩阵的行和列
p=new double[r*c];//为矩阵开辟空间
for(int i=0;i<row*col;i++)//将矩阵元素初始值设为0
p[i]=0;
}
matrix::matrix(int r,int c,double *a)
{
p=new double[r*c];//为矩阵开辟空间
row=r;col=c;//存储矩阵的行和列
for(int i=0;i<row*col;i++)//用数组a初始化矩阵
p[i]=a[i];
}
matrix::matrix(const matrix &a)
{
row=a.row;col=a.col;//存储矩阵的行和列
p=new double[row*col];//为矩阵开辟空间
for(int i=0;i<row*col;i++)//用矩阵a初始化矩阵
p[i]=a.p[i];
}
matrix::~matrix()
{
row=col=0;//更新矩阵的行和列
delete []p;//释放掉原矩阵的空间
p=NULL;//将矩阵指针置空
}
istream& operator>>(istream& in,matrix& x)
{
int r,c;
cout<<"请输入矩阵的行和列："<<endl;
in>>r>>c;
x=matrix(r,c);
cout<<"请输入矩阵的元素值："<<endl;
for(int i=0;i<x.row;i++)
for(int j=0;j<x.col;j++)
in>>x.p[i*x.col+j];
return in;
}
ostream& operator<<(ostream& out,matrix& x)
{
out<<"该矩阵是"<<x.row<<"行"<<x.col<<"列的矩阵"<<endl;
out<<"其元素如下:"<<endl;
for(int i=0;i<x.row;i++)
{
for(int j=0;j<x.col;j++)
{
if(fabs(x.p[i*x.col+j])>=0.0000001)
out<<x.p[i*x.col+j]<<' ';
else cout<<"0 ";
}
out<<endl;
}
return out;
}
matrix &matrix::operator=(const matrix &a)//进行此操作可能会改变原矩阵的规模
{
if(this==&a) return *this;//防止a=a的赋值
row=a.row;//更新矩阵的行
col=a.col;//更新矩阵的列
delete []p;//释放掉原区域
p=new double[row*col];//分配新区域
memcpy(p,a.p,sizeof(double)*row*col);//对矩阵进行赋值
return *this;//返回赋值后的矩阵的引用
}
matrix &matrix::operator+=(const matrix &a)
{
if(row!=a.row||col!=a.col)//矩阵相加要求两个矩阵行数以及列数都要相同
{
cout<<"不满足矩阵相加条件！"<<endl;
return *this;
}
matrix t=(*this)+a;
(*this)=t;
return *this;//返回赋值后的矩阵的引用
}
matrix &matrix::operator-=(const matrix &a)
{
if(row!=a.row||col!=a.col)//矩阵相减要求两个矩阵行数以及列数都要相同
{
cout<<"不满足矩阵相减条件！"<<endl;
return *this;
}
matrix t=(*this)-a;
(*this)=t;
return *this;//返回赋值后的矩阵的引用
}
matrix &matrix::operator*=(const matrix &a)//进行此操作可能会改变原矩阵的规模
{
if(col!=a.row)//矩阵相乘要求左矩阵的列数和右矩阵的行数要相同
{
cout<<"不满足矩阵相乘条件！"<<endl;
return *this;
}
matrix t=(*this)*a;
*this=t;
return *this;//返回赋值后的矩阵的引用
}
matrix matrix::Transpose()
{
matrix t(col,row);
for(int i=0;i<row;i++)//进行元素更新
for(int j=0;j<col;j++)
t.p[j*row+i]=p[i*col+j];
return t;
}
matrix matrix::Transformation()
{
matrix t(*this);
for(int i=0,j=0;i<row&&j<col;i++,j++)//枚举当前处理的行
{
int maxRow=i;//设定当前行为第j列的数的绝对值最大的行
for(int k=i+1;k<row;k++)//寻找当前列绝对值最大的行,减少精度损失
{
if(abs(t.p[k*col+j])>abs(t.p[maxRow*col+j]))
maxRow=k;
}
if(maxRow!=i)//与第row行交换
{
for(int k=i;k<col;k++)
swap(t.p[i*col+k],t.p[maxRow*col+k]);
}
if(fabs(t.p[i*col+j])<0.000001)//col列第row行以下全是0，处理当前行的下一列
{
i--;
continue;
}
for(int k=0;k<row;k++)//枚举要删去的行
{
if(k==i) continue;
if(fabs(t.p[k*col+j])>0.000001)
{
double temp=t.p[k*col+j]/t.p[i*col+j];
for(int l=j;l<col;l++)
t.p[k*col+l]-=t.p[i*col+l]*temp;
}
}
for(int k=j+1;k<col;k++)//将第i行首非零元化为1
t.p[i*col+k]/=t.p[i*col+j];
t.p[i*col+j]=1;
}
return t;
}
double matrix::Det()
{
if(col!=row)//只有方阵才会有行列式
{
cout<<"该矩阵不是方阵，无法计算行列式值!";
return 0;
}
matrix t(*this);
double ans=1;
for(int i=0,j=0;i<row&&j<col;i++,j++)//枚举当前处理的行
{
int maxRow=i;//设定当前行为第j列的数的绝对值最大的行
for(int k=i+1;k<row;k++)//寻找当前列绝对值最大的行,减少精度损失
{
if(abs(t.p[k*col+j])>abs(t.p[maxRow*col+j]))
maxRow=k;
}
if(maxRow!=i)//与第maxRow行交换
{
ans*=-1;//每次交换结果要变号
for(int k=i;k<col;k++)
swap(t.p[i*col+k],t.p[maxRow*col+k]);
}
if(fabs(t.p[i*col+j])<0.000001)//col列第row行以下全是0，处理当前行的下一列
{
i--;
continue;
}
for(int k=i+1;k<row;k++)//枚举要删去的行
{
if(fabs(t.p[k*col+j])>0.000001)
{
double temp=t.p[k*col+j]/t.p[i*col+j];
for(int l=j;l<col;l++)
t.p[k*col+l]-=t.p[i*col+l]*temp;
}
}
}
for(int i=0;i<row;i++)//将ans*对角线上面的元素就得到了行列式的值
ans*=t.p[i*col+i];
return ans;
}
matrix operator+(const matrix &a,const matrix &b)
{
if((a.row!=b.row)||(a.col!=b.col))//矩阵相加要求两个矩阵行数以及列数都要相同
{
cout<<"不满足矩阵相加条件!";
return a;
}
matrix c(a.row,a.col);//利用构造函数初始化一个矩阵c，用来存储答案
int row=a.row;
int col=a.col;
for(int i=0;i<row;i++)//更新矩阵
for(int j=0;j<col;j++)
c.p[i*col+j]=a.p[i*col+j]+b.p[i*col+j];
return c;//返回结果矩阵
}
matrix operator-(const matrix &a,const matrix &b)
{
if((a.row!=b.row)||(a.col!=b.col))//矩阵相减要求两个矩阵行数以及列数都要相同
{
cout<<"不满足矩阵相减条件!";
return a;
}
matrix c(a.row,a.col);//利用构造函数初始化一个矩阵c，用来存储答案
int row=a.row;
int col=a.col;
for(int i=0;i<row;i++)//更新矩阵
for(int j=0;j<col;j++)
c.p[i*col+j]=a.p[i*col+j]-b.p[i*col+j];
return c;//返回结果矩阵
}
matrix operator*(const matrix &a,const matrix &b)
{
if(a.col!=b.row)//矩阵相乘要求左矩阵的列数和右矩阵的行数要相同
{
cout<<"不满足矩阵相乘条件!";
return a;
}
matrix c(a.row,b.col);//结果矩阵的行数与左矩阵的行数相同，列数与右矩阵的列数相同
for(int i=0;i<a.row;i++)//更新矩阵
for(int j=0;j<b.col;j++)
for(int k=0;k<a.col;k++)
c.p[i*c.col+j]+=a.p[i*a.col+k]*b.p[k*b.col+j];
return c;//返回结果矩阵
}
double &matrix::operator()(int r,int c)
{
if(r<0||r>=row||c<0||c>=col)//发生越界
{
cout<<"输入不合法!";
exit(0);
}
return p[r*col+c];//返回结果
}
int matrix::rank()
{
matrix b=Transformation();//先生成一个行初等变换后的矩阵
for(int i=0;i<row;i++)
{
bool flag=true;
for(int j=i;j<col;j++)//如果某行全为0，则在此行之后的行也全为0
{
if(fabs(b.p[i*col+j])>0.000001)//这里是对化成行简化阶梯型的矩阵进行判断
{
flag=false;
break;
}
}
if(flag) return i;//返回非零行信息
}
return row;//满秩矩阵
}
void matrix::cal(matrix b)
{
matrix t(*this);
for(int i=0,j=0;i<row&&j<col;i++,j++)//枚举当前处理的行
{
int maxRow=i;//设定当前行为第j列的数的绝对值最大的行
for(int k=i+1;k<row;k++)//寻找当前列绝对值最大的行,减少精度损失
{
if(abs(t.p[k*col+j])>abs(t.p[maxRow*col+j]))
maxRow=k;
}
if(maxRow!=i)//与第maxRow行交换
{
for(int k=i;k<col;k++)
swap(t.p[i*col+k],t.p[maxRow*col+k]);
swap(b.p[i],b.p[maxRow]);//交换解集
}
if(fabs(t.p[i*col+j])<0.000001)//col列第row行以下全是0，处理当前行的下一列
{
i--;
continue;
}
for(int k=0;k<row;k++)//枚举要删去的行
{
if(k==i) continue;
if(fabs(t.p[k*col+j])>0.000001)
{
double temp=t.p[k*col+j]/t.p[i*col+j];
for(int l=j;l<col;l++)
t.p[k*col+l]-=t.p[i*col+l]*temp;
b.p[k]-=b.p[i]*temp;//更新解集
}
}
}
int cnt=rank();//记录系数矩阵的秩
for(int i=cnt;i<=row;i++)
{
if(fabs(b.p[i])>0.000001)//无解
{
cout<<"方程无解！"<<endl;
return ;
}
}
vector<int> free;//记录自由变元
vector<int> ffree;//记录非自由变元
for(int i=0,j=0;i<cnt;i++,j++)
{
while(fabs(t.p[i*col+j])<0.000001) free.push_back(j++);
ffree.push_back(j);
for(int k=j+1;k<col;k++)
t.p[i*col+k]/=t.p[i*col+j];
b.p[i]/=t.p[i*col+j];
t.p[i*col+j]=1;
}
if(cnt<col)//多解
{
cout<<"方程有多解:"<<endl;
for(int i=ffree[ffree.size()-1]+1;i<col;i++) free.push_back(i);
cout<<"自由元为:";
for(int i=0;i<free.size();i++)
cout<<"x"<<(free[i]+1)<<" ";
cout<<endl<<"解的形式为:"<<endl;
for(int i=0;i<cnt;i++)
{
cout<<"x"<<(ffree[i]+1)<<"=";
cout<<b.p[i];
for(int j=0;j<free.size();j++)
{
if(fabs(t.p[i*col+free[j]])>0.000001)
{
if((-t.p[i*col+free[j]])>0) cout<<"+";
else cout<<"-";
if(fabs(fabs(t.p[i*col+free[j]])-1)>0.000001)
cout<<fabs(t.p[i*col+free[j]]);
cout<<"x"<<(free[j]+1);
}
}
cout<<endl;
}
}
else//唯一解
{
cout<<"方程有唯一解:"<<endl;
for(int i=0;i<col;i++)
{
b.p[i]/=t.p[i*col+i];
cout<<"x"<<(i+1)<<"="<<b.p[i]<<endl;
}
}
return ;
}
matrix matrix::unit(int n)
{
matrix t(n,n);
for(int i=0;i<n;i++)
t.p[i*t.col+i]=1;
return t;
}
matrix matrix::operator^(int n)
{
if(row!=col)
{
cout<<"不满足求幂次条件!"<<endl;
return (*this);
}
matrix t(*this);
matrix E=unit(col);
while(n)
{
if(n&1)
E*=t;
n>>=1;
t*=t ;
}
return E;
}
matrix matrix::aug(const matrix &b)
{
if(row!=b.row)
{
cout<<"无法将两个矩阵合并为增广矩阵!"<<endl;
return (*this);
}
matrix t(row,col+b.col);
for(int i=0;i<row;i++)
for(int j=0;j<col;j++)
t.p[i*(col+b.col)+j]=p[i*col+j];
for(int i=0;i<row;i++)
for(int j=0;j<b.col;j++)
t.p[i*(col+b.col)+col+j]=b.p[i*b.col+j];
return t;
}
matrix matrix::inverse(bool &flag)
{
if(row!=col)
{
flag=false;
cout<<"该矩阵不是方阵，没有逆矩阵!"<<endl;
return (*this);
}
matrix E=unit(row);
matrix t=aug(E);
matrix x=t.Transformation();
int sign=1;
for(int i=0;i<col;i++)
if(fabs(x.p[(row-1)*x.col+i])>=0.000001)
sign=0;
if(sign)
{
flag=false;
cout<<"该矩阵不是满秩矩阵，没有逆矩阵!"<<endl;
return (*this);
}
flag=true;
for(int i=0;i<row;i++)
for(int j=0;j<col;j++)
E.p[i*col+j]=x.p[i*2*col+col+j];
return E;
}
void instruction()
{
cout<<"欢迎进入矩阵类系统\n";
cout<<"***********************************************\n";
cout<<"1：矩阵加法运算\n";
cout<<"2：矩阵减法运算\n";
cout<<"3：矩阵乘法运算\n";
cout<<"4：求矩阵的转置\n";
cout<<"5：求矩阵的行列式\n";
cout<<"6：求矩阵的秩\n";
cout<<"7：求矩阵方程Ax=b\n";
cout<<"8：求矩阵的逆矩阵\n";
cout<<"9:求矩阵的n次幂\n";
cout<<"0：退出\n";
cout<<"***********************************************\n";
}
matrix p1,p2;
int r,c;
int main()
{
int op;
while(1)
{
instruction();
cout<<"请输入您的操作："<<endl;
cin>>op;
switch(op)
{
case 0:
return 0;
case 1:
{
cin>>p1>>p2;
matrix t=p1+p2;
cout<<"执行加法运算后的矩阵："<<endl<<t;
break;
}
case 2:
{
cin>>p1>>p2;
matrix t=p1-p2;
cout<<"执行减法运算后的矩阵："<<endl<<t;
break;
}
case 3:
{
cin>>p1>>p2;
matrix t=p1*p2;
cout<<"执行乘法运算后的矩阵："<<endl<<t;
break;
}
case 4:
{
cin>>p1;
matrix t=p1.Transpose();
cout<<"转置后的矩阵："<<endl<<t;
break;
}
case 5:
{
cin>>p1;
double t=p1.Det();
cout<<"矩阵的行列式的值："<<t<<endl;
break;
}
case 6:
{
cin>>p1;
int t=p1.rank();
cout<<"矩阵的秩为："<<t<<endl;
break;
}
case 7:
{
cin>>p1>>p2;
p1.cal(p2);
break;
}
case 8:
{
cin>>p1;
bool flag;
matrix t=p1.inverse(flag);
if(flag)
cout<<"求逆后的矩阵为："<<endl<<t;
break;
}
case 9:
{
cin>>p1;
cout<<"请输入幂次：";
int n;
cin>>n;
matrix t=p1^n;
cout<<t<<endl;
break;
}
}
system("pause");
system("cls");
}
return 0;
}
扫一扫
340
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
