矩阵迭代法求矩阵的特征值和特征向量

simwe

c语言:
  
void matiteration(pp,mm,xx,zz,nn,ch,itnu,ii,d0)  
//矩阵迭代法求矩阵的特征值和特征向量；  
int nn,*itnu,ii;  
double *pp,*mm,*xx,*zz,*ch,*d0;  
  {double *y,*z1,*a,*b,*mp,*mp1,*d,*d1;  
  int i;  
  z1=calloc(nn,4);  
  d=matmpl(pp,nn,nn,mm,nn);  
  mp=calloc(nn*nn,4);  
  d1=calloc(nn*nn,4);  
  for(i=0;i<nn;i++)  
  {*(z1+i)=*(zz+i);  
  }  
for(i=0;i<nn*nn;i++)  
  {*(d1+i)=*(d0+i);  
  }  
  *itnu=1;  
  if(ii!=0)  
  {a=matmpl(xx,1,nn,mm,nn);  
   mp1=matmpl(a,1,nn,xx,1);  
   b=matmpl(xx,nn,1,xx,nn);  
   mp=matmpl(b,nn,nn,mm,nn);  
   for(i=0;i<nn*nn;i++)  
   {*(mp+i)=*(mp+i)/(*mp1);  
   }  
   d1=matsub(d1,mp,nn);  
   for(i=0;i<nn*nn;i++)  
   {*(d0+i)=*(d1+i);  
   }  
  
   d=matmpl(d,nn,nn,d0,nn);  
  }  
lab1: y=matmpl(d,nn,nn,z1,1);  
  for(i=0;i<nn;i++)  
  {    *(zz+i)=*(y+i)/(*(y+nn-1));  
  }  
  for(i=0;i<nn;i++)  
  { if(fabs(fabs(*(z1+i))-fabs(*(zz+i)))>0.00001)  
  { for(i=0;i<nn;i++)  
  {*(z1+i)=*(zz+i);  
  }  
      *itnu=*itnu+1;  
      goto lab1;  
  }  
  }  
  *ch=(double)fabs(*(y+nn-1));  
  }  
  
  
  
   
double *matmpl(mat1,m1,n1,mat2,n2) //矩阵相乘；  
double *mat1,*mat2;  
int m1,n1,n2;  
{int i,j,k;  
double *mat,c;  
mat=calloc(m1*n2,4);  
for(i=0;i<m1;i++)  
{ for(j=0;j<n2;j++)  
{ c=0;  
   for(k=0;k<n1;k++)  
{ c=*(mat1+i*n1+k)*(*(mat2+k*n2+j))+c;  
}  
    *(mat+i*n2+j)=c;  
  }  
}  
return(mat);  
}  
  
double *matsub(kk,ll,nn) //矩阵相减；  
double *kk,*ll;  
int nn;  
{double *ma;  
int i;  
ma=calloc(nn*nn,4);  
for(i=0;i<nn*nn;i++)  
{*(ma+i)=*(kk+i)-(*(ll+i));  
}  
return(ma);  
}  
  
void matinvert(pp,nn) //矩阵求逆；  
float *pp;  
int nn;  
{float c,d;  
int i,j,m,k,n,ll;  
int *tt;  
ll=0;  
tt=calloc(2*nn,2);  
  n=0;  
for(k=0;k<nn;k++)  
{for(m=k+1;m<nn;m++)  
   { if(m!=(nn-1)&&fabs(*(pp+m*nn+k))>fabs(*(pp+k*nn+k)))  
       {for(j=0;j<nn;j++)  
   { c=*(pp+k*nn+j);  
   *(pp+k*nn+j)=*(pp+m*nn+j);  
   *(pp+m*nn+j)=c;  
   }  
   *(tt+ll)=k;  
   *(tt+ll+1)=m;  
   ll=ll+2;  
       }  
    }  
    d=1/(*(pp+k*nn+k));  
    *(pp+k*nn+k)=d;  
   for(i=0;i<nn;i++)  
   { if(i!=k) *(pp+i*nn+k)=-d*(*(pp+i*nn+k));  
   }  
     for(j=0;j<nn;j++)  
   { if(j!=k)  
   {for(m=0;m<nn;m++)  
   if(m!=k) *(pp+j*nn+m)=*(pp+j*nn+m)+*(pp+j*nn+k)*(*(pp+k*nn+m));  
   }  
   }  
   for(j=0;j<nn;j++)  
   { if(j!=k) *(pp+k*nn+j)=d*(*(pp+k*nn+j));  
   }  
}  
  for(n=ll-2;n>-1;n=n-2)  
{   i=*(tt+n);  
   j=*(tt+n+1);  
   for(m=0;m<nn;m++)  
   {c=*(pp+m*nn+j);  
   *(pp+m*nn+j)=*(pp+m*nn+i);  
   *(pp+m*nn+i)=c;  
   }  
}  
}

fireghg

如果只是用于处理矩阵，我还是觉得用matlab更方便，c的代码可靠性强，但是matlab的数值运算可是更强的。
求特征值：eig();求秩：rank()……
而且c要对使用变量定义，而且对高维数组处理起来并不方便，而且要动态分配地址比较麻烦。
  
不过这个c代码倒是不错。

fireghg

不过，mythworks的matlab就是用c语言写的，所以其实matlab的内核就是c。

aeroplan

    Orr-sormefeld 方程
用maple或者matlab  
解法可以是变分或者谱法
  
先谢了
yanglif@xjtu.edu.cn

yuan

現在看C已經有點暈了﹐