斯密特正交化进行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;  }  }elsepmm=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;  }  }elsepmm=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;  //特征向量最后一行元素置为1midSum = 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");  }