斯密特正交化进行QR分解
矩阵类：
#include<iostream>
#include <stdlib.h>
#include <math.h>
#define iter 60 //迭代次数
using namespace std;
class Matrix
{
public:
Matrix(const Matrix &rhs)//拷贝构造
{
if(this != &rhs) {
row=rhs.row;
col=rhs.col;
size=row*col;
pmm=new double[size];
for(unsigned i=0; i<size; i++)
pmm[i]=rhs.pmm[i];
}
}
Matrix(unsigned r,unsigned c):row(r),col(c)//非方阵
{
size=r*c;
if (size>0)
{
pmm = new double[size];
for(unsigned int j=0;j<size;j++) //init
{
pmm[j]=0.0;
}
}
else
pmm=NULL;
};
Matrix(unsigned n):row(n),col(n)//方阵
{
size=n*n;
if (size>0)
{
pmm = new double[size];
for(unsigned int j=0;j<size;j++) //init
{
pmm[j]=0.0;
}
}
else
pmm=NULL;
};
~Matrix()
{
if (pmm!=NULL)
{
delete []pmm;
}
}
Matrix &operator=(const Matrix& ); //赋值，必须是成员
friend istream &operator>>(istream&, Matrix&);
friend ostream &operator<<(ostream&, Matrix&);
friend Matrix operator+(const Matrix&,const Matrix&);
friend Matrix operator*(const Matrix&,const Matrix&);
double det();
Matrix diag();
void sort();
void VectValue();
unsigned getrow(){return row;}
unsigned getcolumn(){return col;}
double MatrixNorm2();//2范数
double*operator[](int i)//a[i][j]直接访问元素
{
return pmm+i*row;
}
private:
unsigned row,col,size ;
double *pmm;//数组指针
void QR(Matrix&,Matrix&)const;
void EigenVect(const Matrix& eigenValue ,Matrix &eigenVector);
};
//友元仅仅是指定了访问权限，不是一般意义的函数声明，最好在类外再进行一次函数声明。
istream &operator>>(istream &is, Matrix &obj);
ostream &operator<<(ostream &is, Matrix &obj);
//对称性运算符，如算术，相等，应该是普通非成员函数。
Matrix operator*( const Matrix&,const Matrix& );
Matrix operator+( const Matrix&,const Matrix& );
double dets(int n, double *&aa)
{
double *bb = new double[(n-1)*(n-1)];//创建n-1阶的代数余子式阵bb
int p=0,//判断行是否移动
q=0;//代数余子式
if(n == 1)
{
return aa[0];
}
double sum=0.0;//sum为行列式的值
for(int ai=0;ai<n;ai++) // a的行数把矩阵a(nn)赋值到b(n-1)
{
for(int bi=0;bi<n-1;bi++)//把元素aa[ai][0]代数余子式存到bb[][]
{
if(ai>bi) //ai行的代数余子式是：小于ai的行，aa与bb阵，同行赋值
p=0;
else
p=1; //大于等于ai的行，取aa阵的ai+1行赋值给阵bb的bi行
for(int j=0;j<n-1;j++) //从aa的第二列赋值到第n列
{
bb[bi*(n-1)+j]=aa[(bi+p)*n+j+1];
}
}
if(ai%2==0) q=1;//因为列数为0，所以行数是偶数时候，代数余子式为-1.
else q=(-1);
sum+= q* aa[ai*n] *dets(n-1,bb);
}
return sum;
delete []bb;
}
double Matrix::det()
{
if(col==row)
{
return dets(row,pmm) ;
}
else
{
cout<<("行列不相等无法计算")<<endl;
return 0;
}
}
/
istream &operator>>(istream &is, Matrix &obj)
{
for (unsigned i=0; i<obj.size; i++)
{
is>>obj.pmm[i];
}
return is;
}
ostream &operator<<(ostream &out, Matrix &obj)
{
for( unsigned i = 0; i < obj.row; i++) //打印逆矩阵
{
for( unsigned j = 0; j < obj.col; j++)
{
cout<< obj.pmm[i*obj.col+j]<<"\t";
}
cout<<endl;
}
return out;
}
Matrix operator+( const Matrix& lm,const Matrix& rm)
{
Matrix ret(lm.row,lm.col);
for (unsigned i=0;i<ret.size;i++)
{
ret.pmm[i]=lm.pmm[i] + rm.pmm[i];
}
return ret;
}
Matrix operator*( const Matrix& lm,const Matrix& rm)
{
if(lm.size==0||rm.size==0||lm.col != rm.row)
{
Matrix temp(0,0);
temp.pmm=NULL;
return temp; //数据不合法时候，返回空矩阵
}
Matrix ret(lm.row, rm.col);
for( unsigned i=0;i<lm.row;i++)
{
for(unsigned j=0; j< rm.col;j++)
{
for(unsigned k=0; k< lm.col;k++)//lm.col == rm.row
{
ret.pmm[i*rm.col+j]+= lm.pmm[i*lm.col+k] * rm.pmm[k*rm.col+j];
}
}
}
return ret;
}
Matrix& Matrix::operator=( const Matrix& B )
{
row=B.row;
col=B.col;
size=B.size;
for (unsigned i=0;i<size;i++)
{
pmm[i]=B.pmm[i];
}
return *this;
}
//||matrix||_2 求A矩阵的2范数
double Matrix::MatrixNorm2()
{
double norm = 0;
for(unsigned i = 0;i < size;++i)
{
norm += pmm[i] * pmm[i];
}
return (double)sqrt(norm);
}
void Matrix::sort()
{
double tem;
for (unsigned i=0; i<size;i++)
{
for (unsigned j=i+1; j<size;j++)
{
if (pmm[i]<pmm[j])
{
tem=pmm[i];
pmm[i]=pmm[j];
pmm[j]=tem;
}
}
}
}
Matrix Matrix::diag()
{
if (row!=col)
{
Matrix m(0);
return m;
}
Matrix m(row,1);
for (unsigned i=0;i<row;i++)
{
m.pmm[i]=pmm[i*row+i];
}
return m;
}
//----------------------QR分解-------------------------------------------
//将A分解为Q和R
void Matrix::QR(Matrix &Q,Matrix &R) const
{
//如果A不是一个二维方阵，则提示错误，函数计算结束
if(row != col )
{
printf("ERROE: QR() parameter A is not a square matrix!\n");
return;
}
const unsigned N = row;
double *a=new double[N];
double *b=new double[N];
for(unsigned j = 0 ; j < N;++j) //(Gram-Schmidt) 正交化方法
{
for(unsigned i = 0;i < N; ++i) //第j列的数据存到a，b
a[i] = b[i] = pmm[i * N + j];
for(unsigned i = 0; i<j; ++i) //第j列之前的列
{
R.pmm[i * N + j] = 0; //
for(unsigned m = 0;m < N; ++m)
{
R.pmm[i * N + j] += a[m] * Q.pmm[m *N + i]; //R[i,j]值为Q第i列与A的j列的内积
}
for(unsigned m = 0;m < N; ++m)
{
b[m] -= R.pmm[i * N + j] * Q.pmm[m * N + i]; //
}
}
double norm = 0;
for(unsigned i = 0;i < N;++i)
{
norm += b[i] * b[i];
}
norm= (double)sqrt(norm);
R.pmm[j*N + j] = norm; //向量b[]的2范数存到R[j,j]
for(unsigned i = 0; i < N; ++i)
{
Q.pmm[i * N + j] = b[i] / norm; //Q 阵的第j列为单位化的b[]
}
}
delete[]a;
delete[]b;
}
//----------已知矩阵的特征值求矩阵特征向量-----------------
//eigenValue为矩阵A的特征值，
//eigenVector为计算结果即矩阵A的特征向量
void Matrix::EigenVect(const Matrix& eigenValue ,Matrix &eigenVector)
{
const unsigned NUM = col;
double eValue;
double sum,midSum,diag;
Matrix copym(*this);
for(unsigned count = 0; count < NUM;++count)
{
//计算特征值为eValue，求解特征向量时的系数矩阵
*this=copym;
eValue = eigenValue.pmm[count] ;
for(unsigned i = 0 ; i < col ; ++i)//A-lambda*I
{
pmm[i * col + i] -= eValue;
}
//cout<<*this<<endl;
//将 this为阶梯型的上三角矩阵
for(unsigned i = 0 ; i < row - 1 ; ++i)
{
diag = pmm[i*col + i]; //提取对角元素
for(unsigned j = i; j < col; ++j)
{
pmm[i*col + j] /= diag; //【i,i】元素变为1
}
for(unsigned j=i+1; j<row; ++j)
{
diag = pmm[j * col + i];
for(unsigned q = i ; q < col; ++q)//消去第i+1行的第i个元素
{
pmm[j*col + q] -= diag*pmm[i*col + q];
}
}
}
//cout<<*this<<endl;
//特征向量最后一行元素置为1
midSum = eigenVector.pmm[(eigenVector.row - 1) * eigenVector.col + count] = 1;
for(int m = row - 2; m >= 0; --m)
{
sum = 0;
for(unsigned j = m + 1;j < col; ++j)
{
sum += pmm[m * col + j] * eigenVector.pmm[j * eigenVector.col + count];
}
sum= -sum / pmm[m * col + m];
midSum += sum * sum;
eigenVector.pmm[m * eigenVector.col + count] = sum;
}
midSum = sqrt(midSum);
for(unsigned i = 0; i < eigenVector.row ; ++i)
{
eigenVector.pmm[i * eigenVector.col + count] /= midSum; //每次求出一个列向量
}
}
}
void Matrix::VectValue()
{
if (size==0||row!=col)
{
cout<<"矩阵为空或者非方阵！"<<endl;
return;
}
if (det()==0)
{
cout<<"非满秩矩阵没法用QR分解计算特征值！"<<endl;
return;
}
const unsigned N=row;
Matrix A(*this);//备份矩阵
Matrix Q(N),R(N);
/*当迭代次数足够多时,A 趋于上三角矩阵，上三角矩阵的对角元就是A的全部特征值。*/
for(int k = 0;k < iter; ++k)
{
//cout<<"this:\n"<<*this<<endl;
QR(Q,R);
*this=R*Q;
/* cout<<"Q:\n"<<Q<<endl;
cout<<"R:\n"<<R<<endl; */
}
Matrix eigenValue=diag();//特征值
cout<<"\neigenValue:\n"<<eigenValue<<endl;
Matrix eigenVect(N);//特征向量
A.EigenVect (eigenValue,eigenVect);
cout<<"\neigenVect:\n"<<eigenVect<<endl;
}
测试用的主函数：
#include"matrix.h"
#include<iostream>
#include <fstream> // std::ifstream
using namespace std;
int main()
{
cout<<"输入阶数：";
int N=0;
cin>>N;
Matrix mm(N);
int flag=1;
cout<<"屏幕输入矩阵按0,其他数字则从文件mm.txt读入矩阵：";
cin>>flag;
cout<<endl;
if (flag==0)
{
cout<<"input matrix:"<<endl;
cin>>mm;
}else
{
ifstream ifs;
ifs.open ("mm.txt", ifstream::in);
ifs>>mm;
}
mm.VectValue();
system("pause");
} 
