如何用C语言编写求对称矩阵的特征值和特征向量的程序

//数值计算程序-特征值和特征向量////////////////////////////////////////////////////////////////约化对称矩阵为三对角对称矩阵//利用Householder变换将n阶实对称矩阵约化为对称三对角矩阵//a-长度为n*n的数组，存放n阶实对称矩阵//n-矩阵的阶数//q-长度为n*n的数组，返回时存放Householder变换矩阵//b-长度为n的数组，返回时存放三对角阵的主对角线元素//c-长度为n的数组，返回时前n-1个元素存放次对角线元素voideastrq(doublea[],intn,doubleq[],doubleb[],doublec[]);////////////////////////////////////////////////////////////////求实对称三对角对称矩阵的全部特征值及特征向量//利用变型QR方法计算实对称三对角矩阵全部特征值及特征向量//n-矩阵的阶数//b-长度为n的数组，返回时存放三对角阵的主对角线元素//c-长度为n的数组，返回时前n-1个元素存放次对角线元素//q-长度为n*n的数组，若存放单位矩阵，则返回实对称三对角矩阵的特征向量组//若存放Householder变换矩阵，则返回实对称矩阵A的特征向量组//a-长度为n*n的数组，存放n阶实对称矩阵intebstq(intn,doubleb[],doublec[],doubleq[],doubleeps,intl);////////////////////////////////////////////////////////////////约化实矩阵为赫申伯格(Hessenberg)矩阵//利用初等相似变换将n阶实矩阵约化为上H矩阵//a-长度为n*n的数组，存放n阶实矩阵,返回时存放上H矩阵//n-矩阵的阶数voidechbg(doublea[],intn);////////////////////////////////////////////////////////////////求赫申伯格(Hessenberg)矩阵的全部特征值//返回值小于0表示超过迭代jt次仍未达到精度要求//返回值大于0表示正常返回//利用带原点位移的双重步QR方法求上H矩阵的全部特征值//a-长度为n*n的数组，存放上H矩阵//n-矩阵的阶数//u-长度为n的数组，返回n个特征值的实部//v-长度为n的数组，返回n个特征值的虚部//eps-控制精度要求//jt-整型变量，控制最大迭代次数intedqr(doublea[],intn,doubleu[],doublev[],doubleeps,intjt);////////////////////////////////////////////////////////////////求实对称矩阵的特征值及特征向量的雅格比法//利用雅格比(Jacobi)方法求实对称矩阵的全部特征值及特征向量//返回值小于0表示超过迭代jt次仍未达到精度要求//返回值大于0表示正常返回//a-长度为n*n的数组，存放实对称矩阵，返回时对角线存放n个特征值//n-矩阵的阶数//u-长度为n*n的数组，返回特征向量(按列存储)//eps-控制精度要求//jt-整型变量，控制最大迭代次数inteejcb(doublea[],intn,doublev[],doubleeps,intjt);//////////////////////////////////////////////////////////////选自徐世良数值计算程序集(C)每个程序都加上了适当地注释，陆陆续续干了几个月才整理出来的啊。今天都给贴出来了#includestdio.h#includemath.h//约化对称矩阵为三对角对称矩阵//利用Householder变换将n阶实对称矩阵约化为对称三对角矩阵//a-长度为n*n的数组，存放n阶实对称矩阵//n-矩阵的阶数//q-长度为n*n的数组，返回时存放Householder变换矩阵//b-长度为n的数组，返回时存放三对角阵的主对角线元素//c-长度为n的数组，返回时前n-1个元素存放次对角线元素voideastrq(doublea[],intn,doubleq[],doubleb[],doublec[]){inti,j,k,u,v;doubleh,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(i1){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阶实对称矩阵intebstq(intn,doubleb[],doublec[],doubleq[],doubleeps,intl){inti,j,k,m,it,u,v;doubled,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(hd){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);}//约化实矩阵为赫申伯格(Hessenberg)矩阵//利用初等相似变换将n阶实矩阵约化为上H矩阵//a-长度为n*n的数组，存放n阶实矩阵,返回时存放上H矩阵//n-矩阵的阶数voidechbg(doublea[],intn){inti,j,k,u,v;doubled,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;}//求赫申伯格(Hessenberg)矩阵的全部特征值//利用带原点位移的双重步QR方法求上H矩阵的全部特征值//返回值小于0表示超过迭代jt次仍未达到精度要求//返回值大于0表示正常返回//a-长度为n*n的数组，存放上H矩阵//n-矩阵的阶数//u-长度为n的数组，返回n个特征值的实部//v-长度为n的数组，返回n个特征值的虚部//eps-控制精度要求//jt-整型变量，控制最大迭代次数intedqr(doublea[],intn,doubleu[],doublev[],doubleeps,intjt){intm,it,i,j,k,l,ii,jj,kk,ll;doubleb,c,w,g,xy,p,q,r,x,s,e,f,z,y;it=0;m=n;while(m!=0){l=m-1;while((l0)&&(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