1、/数值计算程序-特征值和特征向量 / /约化对称矩阵为三对角对称矩阵 /利用 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); / /求实对称三对角对称矩阵的全部特征值及特征向量 /利用变型
2、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)矩阵 /
3、利用初等相似变换将 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-控制精
4、度要求 /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 eej
5、cb(double a,int n,double v,double eps,int jt); / 选自 每个程序都加上了适当地注释，陆陆续续干了几个月才整理出来的啊。 今天都给贴出来了 #include “stdio.h“ #include “math.h“ /约化对称矩阵为三对角对称矩阵 /利用 Householder 变换将 n 阶实对称矩阵约化为对称三对角矩阵 /a-长度为 n*n 的数组，存放 n 阶实对称矩阵 /n-矩阵的阶数 /q-长度为 n*n 的数组，返回时存放 Householder 变换矩阵 /b-长度为 n 的数组，返回时存放三对角阵的主对角线元素 /c-长度为 n 的数
6、组，返回时前 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=1; i-) h=0.0; if (i1) for (k=0; k0.0) ci-1=-ci-1; h=h-qu*ci-1; qu=qu-ci-1; f=0.0; for (j=0; j=0) for (j=0; j=0) for (j=0; jd) d=h; m=j; while (md) m=m+1; if (m!=j) do if (it=l) p
7、rintf(“failn“); return(-1); it=it+1; g=bj; p=(bj+1-g)/(2.0*cj); r=sqrt(p*p+1.0); if (p=0.0) bj=cj/(p+r); else bj=cj/(p-r); h=g-bj; for (i=j+1; i=j; i-) g=e*ci; h=e*p; if (fabs(p)=fabs(ci) e=ci/p; r=sqrt(e*e+1.0); ci+1=s*p*r; s=e/r; e=1.0/r; else e=p/ci; r=sqrt(e*e+1.0); ci+1=s*ci*r; s=1.0/r; e=e/r;
8、p=e*bi-s*g; bi+1=h+s*(e*g+s*bi); for (k=0; kd); bj=bj+f; for (i=0; ifabs(d) d=t; i=j; if (fabs(d)+1.0!=1.0) if (i!=k) for (j=k-1; j0) 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) um-1=a(m-1)*n+m-1; vm-1=0.0; m=m-1; it=0; else if (l=m-2) b=-(aii+all); c=aii*all-ajj*akk;
9、w=b*b-4.0*c; y=sqrt(fabs(w); if (w0.0) xy=1.0; if (b=jt) printf(“failn“); return(-1); it=it+1; for (j=l+2; j=m-1) j=m-1; for (i=l; ifm) fm=d; p=i; q=j; if (fmjt) return(-1); l=l+1; u=p*n+q; w=p*n+p; t=q*n+p; s=q*n+q; x=-au; y=(as-aw)/2.0; omega=x/sqrt(x*x+y*y); if (y0.0) omega=-omega; sn=1.0+sqrt(1.
10、0-omega*omega); sn=omega/sqrt(2.0*sn); cn=sqrt(1.0-sn*sn); fm=aw; aw=fm*cn*cn+as*sn*sn+au*omega; as=fm*sn*sn+as*cn*cn-au*omega; au=0.0; at=0.0; for (j=0; j=n-1; j+) if (j!=p) w=q*n+j; fm=au; au=fm*cn+aw*sn; aw=-fm*sn+aw*cn; for (i=0; i=n-1; i+) if (i!=p) w=i*n+q; fm=au; au=fm*cn+aw*sn; aw=-fm*sn+aw*cn; for (i=0; i=n-1; i+) u=i*n+p; w=i*n+q; fm=vu; vu=fm*cn+vw*sn; vw=-fm*sn+vw*cn; return(1);
