雅克比求矩阵的特征值特征向量函数改进

来源:互联网 发布:java行业看法 编辑:程序博客网 时间:2024/06/02 17:37
<span style="font-size:18px;">double eejcb(vector<vector<float> > a/*(n)*/,int n/*,vector<vector<double> > v,int jt,int o){     int i,j,p=0,q=0,u,w,t,s,l;    double fm,cn,sn,omega,x,y,d;    l=1;vector <double> b;float sum =0.0;vector<vector<double> > v(n);for(int i=0;i <n;i++) v[i].resize(n);//vector<vector<float> > a(n);for(int i=0;i <n;i++) a[i].resize(n);    for (i=0; i<=n-1; i++){ v[/*i*n+*/i].push_back(1.0);for (j=0; j<=n-1; j++){if (i!=j) {v[i/**n+j*/][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][j]/*[i*n+j]*/);//求绝对值if ((i!=j)&&(d>fm)){ fm=d; p=i; q=j;}}}/*if (fm<eps)  {return l;}*/        if (l>jt)  //?{for (i=0; i<=n-1; i++){// float sum =0.0;sum += fabs(a[i][i]);}for (i=0; i<=n-1; i++){double x = fabs(a[i][i])/sum;b.push_back(x);}return b[o];}        l=l+1;       /* u=p*n+q; w=p*n+p; t=q*n+p; s=q*n+q;*/        x=-a[p][q];y=(a[q][q]-a[p][p])/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[p][p];        a[p][p]=fm*cn*cn+a[q][q]*sn*sn+a[p][q]*omega;        a[q][q]=fm*sn*sn+a[q][q]*cn*cn-a[p][q]*omega;        a[p][q]=0.0;a[q][p]=0.0;        for (j=0; j<=n-1; j++){if ((j!=p)&&(j!=q)){ /*u=p*n+j;w=q*n+j;*/fm=a[p][j];a[p][j]=fm*cn+a[q][j]*sn;a[q][j]=-fm*sn+a[q][j]*cn;}}        for (i=0; i<=n-1; i++){if ((i!=p)&&(i!=q))            { /*u=i*n+p; w=i*n+q;*/fm=a[i][p];a[i][p]=fm*cn+a[i][q]*sn;a[i][q]=-fm*sn+a[i][q]*cn;            }}        for (i=0; i<=n-1; i++){ /* u=i*n+p; w=i*n+q;*/           fm=v[i][p];           v[i][p]=fm*cn+v[i][q]*sn;           v[i][q]=-fm*sn+v[i][q]*cn;}}/*return a[o][o];*///vector <double> b;//float sum =0.0;for (i=0; i<=n-1; i++){// float sum =0.0;sum += fabs(a[i][i]);}for (i=0; i<=n-1; i++){double x = fabs(a[i][i])/sum;b.push_back(x);}return b[o];}</span>
注:网上下载的代码是以一维数组为接口,我将其全部改为二维vector形式
0 0