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; 
   } 
} 
} 

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); 
}