斯密特正交化进行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");  
}