数值计算程序-特征值和特征向量 [转]
2012-04-22 16:58
537 查看
//数值计算程序-特征值和特征向量 ////////////////////////////////////////////////////////////// //约化对称矩阵为三对角对称矩阵 //利用Householder变换将n阶实对称矩阵约化为对称三对角矩阵 //a-长度为n*n的数组,存放n阶实对称矩阵 //n-矩阵的阶数 //q-长度为n*n的数组,返回时存放Householder变换矩阵 //b-长度为n的数组,返回时存放三对角阵的主对角线元素 //c-长度为n的数组,返回时前n-1个元素存放次对角线元素 void eastrq(double a[],int n,double q[],double b[],double c[]); ////////////////////////////////////////////////////////////// //求实对称三对角对称矩阵的全部特征值及特征向量 //利用变型QR方法计算实对称三对角矩阵全部特征值及特征向量 //n-矩阵的阶数 //b-长度为n的数组,返回时存放三对角阵的主对角线元素 //c-长度为n的数组,返回时前n-1个元素存放次对角线元素 //q-长度为n*n的数组,若存放单位矩阵,则返回实对称三对角矩阵的特征向量组 // 若存放Householder变换矩阵,则返回实对称矩阵A的特征向量组 //a-长度为n*n的数组,存放n阶实对称矩阵 int ebstq(int n,double b[],double c[],double q[],double eps,int l); ////////////////////////////////////////////////////////////// //约化实矩阵为赫申伯格(Hessen berg)矩阵 //利用初等相似变换将n阶实矩阵约化为上H矩阵 //a-长度为n*n的数组,存放n阶实矩阵,返回时存放上H矩阵 //n-矩阵的阶数 void echbg(double a[],int n); ////////////////////////////////////////////////////////////// //求赫申伯格(Hessen berg)矩阵的全部特征值 //返回值小于0表示超过迭代jt次仍未达到精度要求 //返回值大于0表示正常返回 //利用带原点位移的双重步QR方法求上H矩阵的全部特征值 //a-长度为n*n的数组,存放上H矩阵 //n-矩阵的阶数 //u-长度为n的数组,返回n个特征值的实部 //v-长度为n的数组,返回n个特征值的虚部 //eps-控制精度要求 //jt-整型变量,控制最大迭代次数 int edqr(double a[],int n,double u[],double v[],double eps,int jt); ////////////////////////////////////////////////////////////// //求实对称矩阵的特征值及特征向量的雅格比法 //利用雅格比(Jacobi)方法求实对称矩阵的全部特征值及特征向量 //返回值小于0表示超过迭代jt次仍未达到精度要求 //返回值大于0表示正常返回 //a-长度为n*n的数组,存放实对称矩阵,返回时对角线存放n个特征值 //n-矩阵的阶数 //u-长度为n*n的数组,返回特征向量(按列存储) //eps-控制精度要求 //jt-整型变量,控制最大迭代次数 int eejcb(double a[],int n,double v[],double eps,int jt); ////////////////////////////////////////////////////////////// 选自<<徐世良数值计算程序集(C)>> 每个程序都加上了适当地注释,陆陆续续干了几个月才整理出来的啊。 今天都给贴出来了 #include "stdio.h" #include "math.h" //约化对称矩阵为三对角对称矩阵 //利用Householder变换将n阶实对称矩阵约化为对称三对角矩阵 //a-长度为n*n的数组,存放n阶实对称矩阵 //n-矩阵的阶数 //q-长度为n*n的数组,返回时存放Householder变换矩阵 //b-长度为n的数组,返回时存放三对角阵的主对角线元素 //c-长度为n的数组,返回时前n-1个元素存放次对角线元素 void eastrq(double a[],int n,double q[],double b[],double c[]) { int i,j,k,u,v; double h,f,g,h2; for (i=0; i<=n-1; i++) { for (j=0; j<=n-1; j++) { u=i*n+j; q[u]=a[u]; } } for (i=n-1; i>=1; i--) { h=0.0; if (i>1) { for (k=0; k<=i-1; k++) { u=i*n+k; h=h+q[u]*q[u]; } } if (h+1.0==1.0) { c[i-1]=0.0; if (i==1) { c[i-1]=q[i*n+i-1]; } b[i]=0.0; } else { c[i-1]=sqrt(h); u=i*n+i-1; if (q[u]>0.0) { c[i-1]=-c[i-1]; } h=h-q[u]*c[i-1]; q[u]=q[u]-c[i-1]; f=0.0; for (j=0; j<=i-1; j++) { q[j*n+i]=q[i*n+j]/h; g=0.0; for (k=0; k<=j; k++) { g=g+q[j*n+k]*q[i*n+k]; } if (j+1<=i-1) { for (k=j+1; k<=i-1; k++) { g=g+q[k*n+j]*q[i*n+k]; } } c[j-1]=g/h; f=f+g*q[j*n+i]; } h2=f/(h+h); for (j=0; j<=i-1; j++) { f=q[i*n+j]; g=c[j-1]-h2*f; c[j-1]=g; for (k=0; k<=j; k++) { u=j*n+k; q[u]=q[u]-f*c[k-1]-g*q[i*n+k]; } } b[i]=h; } } b[0]=0.0; for (i=0; i<=n-1; i++) { if ((b[i]!=0.0)&&(i-1>=0)) { for (j=0; j<=i-1; j++) { g=0.0; for (k=0; k<=i-1; k++) { g=g+q[i*n+k]*q[k*n+j]; } for (k=0; k<=i-1; k++) { u=k*n+j; q[u]=q[u]-g*q[k*n+i]; } } } u=i*n+i; b[i]=q[u]; q[u]=1.0; if (i-1>=0) { for (j=0; j<=i-1; j++) { q[i*n+j]=0.0; q[j*n+i]=0.0; } } } return; //求实对称三对角对称矩阵的全部特征值及特征向量 //利用变型QR方法计算实对称三对角矩阵全部特征值及特征向量 //n-矩阵的阶数 //b-长度为n的数组,返回时存放三对角阵的主对角线元素 //c-长度为n的数组,返回时前n-1个元素存放次对角线元素 //q-长度为n*n的数组,若存放单位矩阵,则返回实对称三对角矩阵的特征向量组 // 若存放Householder变换矩阵,则返回实对称矩阵A的特征向量组 //a-长度为n*n的数组,存放n阶实对称矩阵 int ebstq(int n,double b[],double c[],double q[],double eps,int l) { int i,j,k,m,it,u,v; double d,f,h,g,p,r,e,s; c[n-1]=0.0; d=0.0; f=0.0; for (j=0; j<=n-1; j++) { it=0; h=eps*(fabs(b[j])+fabs(c[j])); if (h>d) { d=h; } m=j; while ((m<=n-1)&&(fabs(c[m])>d)) { m=m+1; } if (m!=j) { do { if (it==l) { printf("fail\n"); return(-1); } it=it+1; g=b[j]; p=(b[j+1]-g)/(2.0*c[j]); r=sqrt(p*p+1.0); if (p>=0.0) { b[j]=c[j]/(p+r); } else { b[j]=c[j]/(p-r); } h=g-b[j]; for (i=j+1; i<=n-1; i++) { b[i]=b[i]-h; } f=f+h; p=b[m]; e=1.0; s=0.0; for (i=m-1; i>=j; i--) { g=e*c[i]; h=e*p; if (fabs(p)>=fabs(c[i])) { e=c[i]/p; r=sqrt(e*e+1.0); c[i+1]=s*p*r; s=e/r; e=1.0/r; } else { e=p/c[i]; r=sqrt(e*e+1.0); c[i+1]=s*c[i]*r; s=1.0/r; e=e/r; } p=e*b[i]-s*g; b[i+1]=h+s*(e*g+s*b[i]); for (k=0; k<=n-1; k++) { u=k*n+i+1; v=u-1; h=q[u]; q[u]=s*q[v]+e*h; q[v]=e*q[v]-s*h; } } c[j]=s*p; b[j]=e*p; } while (fabs(c[j])>d); } b[j]=b[j]+f; } for (i=0; i<=n-1; i++) { k=i; p=b[i]; if (i+1<=n-1) { j=i+1; while ((j<=n-1)&&(b[j]<=p)) { k=j; p=b[j]; j=j+1; } } if (k!=i) { b[k]=b[i]; b[i]=p; for (j=0; j<=n-1; j++) { u=j*n+i; v=j*n+k; p=q[u]; q[u]=q[v]; q[v]=p; } } } return(1); } //约化实矩阵为赫申伯格(Hessen berg)矩阵 //利用初等相似变换将n阶实矩阵约化为上H矩阵 //a-长度为n*n的数组,存放n阶实矩阵,返回时存放上H矩阵 //n-矩阵的阶数 void echbg(double a[],int n) { int i,j,k,u,v; double d,t; for (k=1; k<=n-2; k++) { d=0.0; for (j=k; j<=n-1; j++) { u=j*n+k-1; t=a[u]; if (fabs(t)>fabs(d)) { d=t; i=j; } } if (fabs(d)+1.0!=1.0) { if (i!=k) { for (j=k-1; j<=n-1; j++) { u=i*n+j; v=k*n+j; t=a[u]; a[u]=a[v]; a[v]=t; } for (j=0; j<=n-1; j++) { u=j*n+i; v=j*n+k; t=a[u]; a[u]=a[v]; a[v]=t; } } for (i=k+1; i<=n-1; i++) { u=i*n+k-1; t=a[u]/d; a[u]=0.0; for (j=k; j<=n-1; j++) { v=i*n+j; a[v]=a[v]-t*a[k*n+j]; } for (j=0; j<=n-1; j++) { v=j*n+k; a[v]=a[v]+t*a[j*n+i]; } } } } return; } //求赫申伯格(Hessen berg)矩阵的全部特征值 //利用带原点位移的双重步QR方法求上H矩阵的全部特征值 //返回值小于0表示超过迭代jt次仍未达到精度要求 //返回值大于0表示正常返回 //a-长度为n*n的数组,存放上H矩阵 //n-矩阵的阶数 //u-长度为n的数组,返回n个特征值的实部 //v-长度为n的数组,返回n个特征值的虚部 //eps-控制精度要求 //jt-整型变量,控制最大迭代次数 int edqr(double a[],int n,double u[],double v[],double eps,int jt) { int m,it,i,j,k,l,ii,jj,kk,ll; double b,c,w,g,xy,p,q,r,x,s,e,f,z,y; it=0; m=n; while (m!=0) { l=m-1; while ((l>0)&&(fabs(a[l*n+l-1])>eps*(fabs(a[(l-1)*n+l-1])+fabs(a[l*n+l])))) { l=l-1; } ii=(m-1)*n+m-1; jj=(m-1)*n+m-2; kk=(m-2)*n+m-1; ll=(m-2)*n+m-2; if (l==m-1) { u[m-1]=a[(m-1)*n+m-1]; v[m-1]=0.0; m=m-1; it=0; } else if (l==m-2) { b=-(a[ii]+a[ll]); c=a[ii]*a[ll]-a[jj]*a[kk]; w=b*b-4.0*c; y=sqrt(fabs(w)); if (w>0.0) { xy=1.0; if (b<0.0) { xy=-1.0; } u[m-1]=(-b-xy*y)/2.0; u[m-2]=c/u[m-1]; v[m-1]=0.0; v[m-2]=0.0; } else { u[m-1]=-b/2.0; u[m-2]=u[m-1]; v[m-1]=y/2.0; v[m-2]=-v[m-1]; } m=m-2; it=0; } else { if (it>=jt) { printf("fail\n"); return(-1); } it=it+1; for (j=l+2; j<=m-1; j++) { a[j*n+j-2]=0.0; } for (j=l+3; j<=m-1; j++) { a[j*n+j-3]=0.0; } for (k=l; k<=m-2; k++) { if (k!=l) { p=a[k*n+k-1]; q=a[(k+1)*n+k-1]; r=0.0; if (k!=m-2) { r=a[(k+2)*n+k-1]; } } else { x=a[ii]+a[ll]; y=a[ll]*a[ii]-a[kk]*a[jj]; ii=l*n+l; jj=l*n+l+1; kk=(l+1)*n+l; ll=(l+1)*n+l+1; p=a[ii]*(a[ii]-x)+a[jj]*a[kk]+y; q=a[kk]*(a[ii]+a[ll]-x); r=a[kk]*a[(l+2)*n+l+1]; } if ((fabs(p)+fabs(q)+fabs(r))!=0.0) { xy=1.0; if (p<0.0) { xy=-1.0; } s=xy*sqrt(p*p+q*q+r*r); if (k!=l) { a[k*n+k-1]=-s; } e=-q/s; f=-r/s; x=-p/s; y=-x-f*r/(p+s); g=e*r/(p+s); z=-x-e*q/(p+s); for (j=k; j<=m-1; j++) { ii=k*n+j; jj=(k+1)*n+j; p=x*a[ii]+e*a[jj]; q=e*a[ii]+y*a[jj]; r=f*a[ii]+g*a[jj]; if (k!=m-2) { kk=(k+2)*n+j; p=p+f*a[kk]; q=q+g*a[kk]; r=r+z*a[kk]; a[kk]=r; } a[jj]=q; a[ii]=p; } j=k+3; if (j>=m-1) { j=m-1; } for (i=l; i<=j; i++) { ii=i*n+k; jj=i*n+k+1; p=x*a[ii]+e*a[jj]; q=e*a[ii]+y*a[jj]; r=f*a[ii]+g*a[jj]; if (k!=m-2) { kk=i*n+k+2; p=p+f*a[kk]; q=q+g*a[kk]; r=r+z*a[kk]; a[kk]=r; } a[jj]=q; a[ii]=p; } } } } } return(1); } //求实对称矩阵的特征值及特征向量的雅格比法 //利用雅格比(Jacobi)方法求实对称矩阵的全部特征值及特征向量 //返回值小于0表示超过迭代jt次仍未达到精度要求 //返回值大于0表示正常返回 //a-长度为n*n的数组,存放实对称矩阵,返回时对角线存放n个特征值 //n-矩阵的阶数 //u-长度为n*n的数组,返回特征向量(按列存储) //eps-控制精度要求 //jt-整型变量,控制最大迭代次数 int eejcb(double a[],int n,double v[],double eps,int jt) { int i,j,p,q,u,w,t,s,l; double fm,cn,sn,omega,x,y,d; l=1; for (i=0; i<=n-1; i++) { v[i*n+i]=1.0; for (j=0; j<=n-1; j++) { if (i!=j) { v[i*n+j]=0.0; } } } while (1==1) { fm=0.0; for (i=0; i<=n-1; i++) { for (j=0; j<=n-1; j++) { d=fabs(a[i*n+j]); if ((i!=j)&&(d>fm)) { fm=d; p=i; q=j; } } } if (fm<eps) { return(1); } if (l>jt) { return(-1); } l=l+1; u=p*n+q; w=p*n+p; t=q*n+p; s=q*n+q; x=-a[u]; y=(a[s]-a[w])/2.0; omega=x/sqrt(x*x+y*y); if (y<0.0) { omega=-omega; } sn=1.0+sqrt(1.0-omega*omega); sn=omega/sqrt(2.0*sn); cn=sqrt(1.0-sn*sn); fm=a[w]; a[w]=fm*cn*cn+a[s]*sn*sn+a[u]*omega; a[s]=fm*sn*sn+a[s]*cn*cn-a[u]*omega; a[u]=0.0; a[t]=0.0; for (j=0; j<=n-1; j++) { if ((j!=p)&&(j!=q)) { u=p*n+j; w=q*n+j; fm=a[u]; a[u]=fm*cn+a[w]*sn; a[w]=-fm*sn+a[w]*cn; } } for (i=0; i<=n-1; i++) { if ((i!=p)&&(i!=q)) { u=i*n+p; w=i*n+q; fm=a[u]; a[u]=fm*cn+a[w]*sn; a[w]=-fm*sn+a[w]*cn; } } for (i=0; i<=n-1; i++) { u=i*n+p; w=i*n+q; fm=v[u]; v[u]=fm*cn+v[w]*sn; v[w]=-fm*sn+v[w]*cn; } } return(1); }
相关文章推荐
- 数值计算程序-特征值和特征向量 [转]
- 特征值与特征向量的数值计算; Matrix Eigenvalues and Eigenvectors Calculating
- 数值分析 幂法求矩阵A按模最大的特征值和相应的特征向量
- C语言调用MKL进行矩阵计算—特征值特征向量求解
- 线性代数(二十九) :特征值与特征向量的计算
- Jacobi 方法计算特征值特征向量
- 理解不同特征值的特征向量之间的非线性关系
- 个人喜欢的关于模式识别、机器学习、推荐系统、图像特征、数值计算、目标跟踪等方面个人主页及博客
- 特征值与特征向量
- MIT18.06线性代数课程笔记21:特征向量、特征值简介、求法与性质
- 矩阵的特征值与特征向量以及相似矩阵
- 均值、方差、协方差、协方差矩阵、特征值、特征向量
- 再谈特征值与特征向量的求取
- 实对称阵属于不同特征值的的特征向量是正交的
- CUDA计算向量内积的程序(源自CUDA范例编程)
- 幂法求矩阵的最大特征值和对应特征向量
- 特征值与特征向量
- Haar-like矩形特征的特征值的快速计算方法
- 矩阵及其变换、特征值与特征向量的物理意义
- 求特征值 特征向量1