QRR.rar_qr分解_qr分解C_矩阵正交

//Eigen_QR.h Householder+QR求解实对称矩阵全部特征值和特征向量 //作者 付阳 南昌大学 //最后修改时间 2009,7 //Email:fyten@163.com #ifndef EIGEN_HQR_H #define EIGEN_HQR_H #include <math.h> #include <omp.h> template <typename Real> class Eigen_QR { private: Real* matrix; Real* work; Real* evalue; int n; Real eps; public: Eigen_QR(Real* cov,Real*gwork,Real* gvalue,int nDim,Real ieps); ~Eigen_QR(void); void EigenHouseholderQR(); void EigenHouseholder(); void EigenHouseholderTransformation(); void HouseholderTransformation(Real h,const int i,const int nn); void HouseholderTransformationNotEqual(Real h,const int i,const int nn); void TransitionCompute(const Real h,const int i,const int nn); Real Computef(const Real h1,const int i,const int nn); void ComputeTransition(const int i,const Real h2,const int nn); void MatrixTridiagonal(); void ZeroSet(const int i); void EigenQR(); bool QRTransformation(const int m,const int j,const Real d,Real& valueoffset); void ComputeValueoffset(const int j,Real& valueoffset); void MatrixElementCompute(const int m,const int j); void ComputeEandS(const Real p,const int i,Real& e,Real& s); void ComputeEandSAbove(const Real p,const int i,Real& e,Real& s); void ComputeEandSBelow(const Real p,const int i,Real& e,Real& s); void EigenVectorCompute(const int i,const Real e,const Real s); void MatrixSwap(const int i,const int j); int Quick(const int u,const int v); void QuickSort(const int u,const int v); void GetValue(Real* gvalue); Real gSum(const int i,const int j,const int nn); Real ComputeSumAbove(const int i,const int j,const int nn,const int jn); Real ComputeSumBelow(const int i,const int j,const int nn); void MatrixIJCompute(const int i,const int j,const Real g); void GetMatrix(Real* gmatrix); void QRLoop(); void QRLoopBar(const int j,Real& d); void QRLoopTransformation(const int j,const Real d,Real& valueoffset); void GetWork(Real* gwork); void MatrixElementSwap(); inline bool FloatEqual(float lhs, float rhs) { if (fabs(lhs - rhs) >= 1.0e-6F) return false; else return true; } inline bool FloatEqual(double lhs, double rhs) { if (fabs(lhs - rhs) >= 1.0e-15) return false; else return true; } inline bool FloatEqual(long double lhs, long double rhs) { if (fabs(lhs - rhs) >= 1.0e-30) return false; else return true; } }; template <typename Real> Eigen_QR<Real>::Eigen_QR(Real* cov,Real*gwork,Real* gvalue,int nDim,Real ieps):matrix(cov),work(gwork),evalue(gvalue),n(nDim),eps(ieps) { } template <typename Real> Eigen_QR<Real>::~Eigen_QR(void) { } template <typename Real> void Eigen_QR<Real>::EigenHouseholderQR() { EigenHouseholder(); EigenQR(); } template <typename Real> void Eigen_QR<Real>::EigenHouseholder() { EigenHouseholderTransformation(); for(int i=0; i<n-1; i++) work[i]=work[i+1]; work[n-1]=0.0; evalue[0]=0.0; MatrixTridiagonal(); MatrixElementSwap(); } template <typename Real> void Eigen_QR<Real>::EigenHouseholderTransformation() { for(int i=n-1; i>=1; i--) { Real h=0.0; int nn = i*n; if(i>1) for(int k=0; k<i; k++) { h += matrix[nn+k]*matrix[nn+k]; } HouseholderTransformation(h,i,nn); } } template <typename Real> void Eigen_QR<Real>::HouseholderTransformation(Real h,const int i,const int nn) { if(FloatEqual(h,0)) { work[i]=0.0; if(i==1) work[i]=matrix[nn+i-1]; evalue[i]=0.0; } else { HouseholderTransformationNotEqual(h,i,nn); } } template <typename Real> void Eigen_QR<Real>::HouseholderTransformationNotEqual(Real h,const int i,const int nn) { work[i]=sqrt(h); if(matrix[nn+i-1]>0.0) work[i]=-work[i]; h=h-matrix[nn+i-1]*work[i]; matrix[nn+i-1]=matrix[nn+i-1]-work[i]; TransitionCompute(h,i,nn); evalue[i]=h; } template <typename Real> void Eigen_QR<Real>::TransitionCompute(const Real h,const int i,const int nn) { Real h1 = 1 / h; Real f=Computef(h1,i,nn); Real h2 = f / (h+h); ComputeTransition(i,h2,nn); } template <typename Real> Real Eigen_QR<Real>::Computef(const Real h1,const int i,const int nn) { Real f=0.0; for(int j=0; j<i; j++) { int jn = j*n; matrix[jn+i]= matrix[nn+j] * h1; Real g = ComputeSumAbove(i,j,nn,jn) + ComputeSumBelow(i,j,nn); work[j] = g * h1; f = f + g * matrix[jn+i]; } return f; } template <typename Real> Real Eigen_QR<Real>::ComputeSumAbove(const int i,const int j,const int nn,const int jn) { Real g = 0.0; for(int k=0; k<=j; k++) { g += matrix[jn+k] * matrix[nn+k]; } return g; } template <typename Real> Real Eigen_QR<Real>::ComputeSumBelow(const int i,const int j,const int nn) { Real g = 0.0; for(int k=j+1; k<i; k++) { int kn = k*n; g += matrix[kn+j] * matrix[nn+k]; } return g; } template <typename Real> void Eigen_QR<Real>::ComputeTransition(const int i,const Real h2,const int nn) { for(int j=0; j<i; j++) { Real f = matrix[nn+j]; Real g = work[j] -h2 * f; work[j] = g; int jn = j*n; for(int k=0; k<=j; k++) matrix[jn+k]=matrix[jn+k]-f*work[k]-g*matrix[nn+k]; } } template <typename Real> void Eigen_QR<Real>::MatrixTridiagonal() { for(int i=0; i<n; i++) { int nn = i*n; if((evalue[i]!=0.0)&&(i-1>=0)) #pragma omp parallel for for(int j=0; j<i; j++) { Real g=gSum(i,j,nn); MatrixIJCompute(i,j,g); } evalue[i]=matrix[nn+i]; matrix[nn+i]=1.0; if(i-1>=0) ZeroSet(i); } } template <typename Real> Real Eigen_QR<Real>::gSum(const int i,const int j,const int nn) { Real g = 0.0; for(int x=0; x<i; x++) { int xn = x*n; g+=matrix[nn+x]*matrix[xn+j]; } return g; } template <typename Real> void Eigen_QR<Real>::MatrixIJCompute(const int i,const int j,const Real g) { int nn = i*n; for(Real *m1 = matrix; m1<matrix+nn; m1+=n) { m1[j] -= g*m1[i]; } } template <typename Real> void Eigen_QR<Real>::ZeroSet(const int i) { for(int j=0; j<i; j++) { matrix[i*n+j]=0.0; matrix[j*n+i]=0.0; } } template <typename Real> void Eigen_QR<Real>::EigenQR() { QRLoop(); QuickSort(0,n-1); } template <typename Real> void Eigen_QR<Real>::QRLoop() { Real d = 0.0; work[n-1]=0.0; Real valueoffset = 0.0; for(int j=0; j<n; j++)//一个j算出一个特征值和其对应的特征向量 { QRLoopBar(j,d); QRLoopTransformation(j,d,valueoffset); } } template <typename Real> void Eigen_QR<Real>::QRLoopTransformation(const int j,const Real d,Real& valueoffset) { int m=j; while((m<n)&&(fabs(work[m])>d)) m++; if(m!=j) { while(QRTransformation(m,j,d,valueoffset)); } evalue[j]+=valueoffset; } template <typename Real> void Eigen_QR<Real>::QRLoopBar(const int j,Real& d) { Real h=eps*(fabs(evalue[j])+fabs(work[j])); if(h>d) d=h; } template <typename Real> bool Eigen_QR<Real>::QRTransformation(const int m,const int j,const Real d,Real& valueoffset) { ComputeValueoffset(j,valueoffset); MatrixElementCompute(m,j); return fabs(work[j])>d; } template <typename Real> void Eigen_QR<Real>::ComputeValueoffset(const int j,Real& valueoffset) { Real g=evalue[j]; Real p=(evalue[j+1]-g)/(2.0*work[j]); Real r=_hypot(p,1.0); if (p<0) { r = -r; } evalue[j]=work[j]/(p+r); Real h=g-evalue[j]; for(int i=j+1; i<n; i++) { evalue[i]-=h; } valueoffset+=h; //valueoffset以后没有变化了 } template <typename Real> void Eigen_QR<Real>::MatrixElementCompute(const int m,const int j) { Real p=evalue[m]; Real e = 1.0; Real s = 0.0; for(int i=m-1; i>=j; i--) { Real g=e*work[i];//work 依然是householder算法传下来的work Real h=e*p; ComputeEandS(p,i,e,s); p=e*evalue[i]-s*g; evalue[i+1]=h+s*(e*g+s*evalue[i]); EigenVectorCompute(i,e,s); } work[j]=s*p; evalue[j]=e*p; } template <typename Real> void Eigen_QR<Real>::ComputeEandS(