fanmifa.cpp

反幂法求矩阵的特征值和特征向量 反幂法适用于矩阵的按模最小的特征值和对应的特征向量。 使用说明： 一般的使用过程： 1、修改输入数据 input2.txt 2、编辑源文件

#include <stdio.h>
#include <stdlib.h>
#include <conio.h>
#include <time.h>
#include <math.h>
#ifndef EPSILON
#define EPSILON 1.0e-12
#endif
double invpower2(int init, double *, double * ,int);
void power( int fs, double *, double *,int );
double lamb1;
double *V;

int symeigen(double eigen[],double A[],double V0[],int N)
{
int i,l,j,step,rank;
static double mlamb;
double *B;

double *x0;
B=(double *)calloc(N*N,sizeof(double));
V=(double *)calloc(N*N,sizeof(double));
x0=(double *)calloc(N,sizeof(double));
for(i=0;i<N;i++)
  { eigen[i]=0.0; for(j=0;j<N;j++) V[i*N+j]=0.0; }
rank=0;
for(step=0;step<N;step++)
  {
  power(step,x0,A,N);
  if(lamb1==0.0) { eigen[step]=0.0; goto bbb; }
  mlamb=lamb1;
  for(i=0;i<N;i++)
    {
    for(j=0;j<N;j++) B[i*N+j]=A[i*N+j];
    B[i*N+i]=A[i*N+i]-lamb1;
    }
  lamb1=invpower2( 1, x0, B ,N);
  if( fabs(lamb1) > EPSILON ) eigen[step]=mlamb+1.0/lamb1;
  if(lamb1==0.0) eigen[step]=mlamb;
bbb:if( fabs(eigen[step]) > EPSILON ) rank=rank+1;
  for(i=0;i<N;i++) V[i*N+step]=x0[i];
  }
for(i=0;i<N;i++)
   for(j=0;j<N;j++) V0[i*N+j]=V[i*N+j];
free(V);
free(B);
free(x0);
return(rank);
}

void power(int i0,double x0[],double A[],int N)
{
static double lamb,sum,err,err1,eps;
int i,j,jj,ii;
double *x1;
time_t t;
x1=(double *)calloc(N*N,sizeof(double));
eps=EPSILON*1000.0;
raa:ii=0;
srand(t);                                         /*for creat x0[]*/
for(j=0;j<i0;j++) { for(i=0;i<N;i++) x0[0]=rand(); }
sum=0.0;
for(i=0;i<N;i++)
  { x0[i]=rand(); sum=sum+x0[i]*x0[i]; }
sum=sqrt(sum);
if( sum < EPSILON ) goto raa;
for(i=0;i<N;i++) x1[i]=x0[i]/sum;      /* x0[] is created*/
lamb=0.5;
goto bbb;
aaa:ii=ii+1;
    for(j=0;j<N;j++)
      {
      x1[j]=0.0;
      for(jj=0;jj<N;jj++) x1[j]=x1[j]+A[j*N+jj]*x0[jj];
      }
    lamb1=0.0;
    for(j=0;j<N;j++) lamb1=lamb1+x0[j]*x1[j];
    if( fabs(lamb1) < EPSILON )
       {
       lamb1=0.0; err=0.0; err1=0.0;
       for(j=0;j<N;j++) x1[j]=x0[j];
       goto bbb;
       }
    else
       { err=fabs(lamb1-lamb); err1=err/fabs(lamb1); }
    lamb=lamb1;
bbb:if(i0>0)
      {
      for(i=0;i<i0;i++)
        {
        sum=0.0;
        for(j=0;j<N;j++) sum=sum+x1[j]*V[j*N+i];
        for(j=0;j<N;j++) x1[j]=x1[j]-sum*V[j*N+i];
        }
      }
sum=0.0;
for(j=0;j<N;j++) sum=sum+x1[j]*x1[j];
sum=sqrt(sum);
if( sum < EPSILON ) goto endd;
for(j=0;j<N;j++) x0[j]=x1[j]/sum;
if( ii < 10 ) goto aaa;
if( ( err > eps ) || (err1 > eps) ) goto aaa;
free(x1);
endd:;
}

/* invpower2.c: Inverse power method for eigenvalues & eigenvectors. */
double invpower2( int init, double x0[], double B[],int N )
{
int i,j,k,i0,istep;
static double lamb,lamb1,sum,err,temp;
int *ik;
double *x1;
time_t t;
ik=(int *)calloc(N,sizeof(int));
x1=(double *)calloc(N*N,sizeof(double));
istep=0;
if( init == 0 )              /* for creat x0[] */
  {
  srand(t);
raa: sum=0.0;
  for(i=0;i<N;i++)
    { x0[i]=rand(); sum=sum+x0[i]*x0[i]; }
  sum=sqrt(sum);
  if( sum < EPSILON ) goto raa;
  for(i=0;i<N;i++) x0[i]=x0[i]/sum;      /* x0[] is created*/
  }
for(k=0;k<N;k++)                 /*colum pivot Doolittle decompsition*/
  {
  temp=fabs( B[k*N+k] );       ik[k]=k;
  for(i=k;i<N;i++)
    if( fabs(B[i*N+k]) > temp ) { temp=fabs(B[i*N+k]); ik[k]=i; }
  i0=ik[k];
  if( i0 != k )
    {
    for(j=0;j<N;j++)
      { temp=B[k*N+j]; B[k*N+j]=B[i0*N+j]; B[i0*N+j]=temp; }
    }
  if( fabs(B[k*N+k]) < EPSILON ) { lamb1=0.0; goto end; }
  for(i=k+1;i<N;i++)
    {
    B[i*N+k]=B[i*N+k]/B[k*N+k];
    for(j=k+1;j<N;j++) B[i*N+j]=B[i*N+j]-B[i*N+k]*B[k*N+j];
    }
  }           /* decompsition is ended */
aaa:istep=istep+1;
    lamb=lamb1;
    for(j=0;j<N;j++) x1[j]=x0[j];
    for(i=0;i<N;i++)                             /* solut LR=Px0[] */
      { j=ik[i]; temp=x0[i]; x0[i]=x0[j]; x0[j]=temp; }
    for(i=1;i<N;i++)
      for(j=0;j<i;j++) x0[i]=x0[i]-B[i*N+j]*x0[j];
    for(i=N-1;i>=0;i--)
      {
      for(j=i+1;j<N;j++) x0[i]=x0[i]-B[i*N+j]*x0[j];
      x0[i]=x0[i]/B[i*N+i];
      }                                         /* end for solution */
    lamb1=0.0;
    for(j=0;j<N;j++) lamb1=lamb1+x0[j]*x1[j];
    err=fabs(lamb1-lamb)/(fabs(lamb1)+1.0);
    sum=0.0;
    for(j=0;j<N;j++) sum=sum+x0[j]*x0[j];
    sum=sqrt(sum);
    for(j=0;j<N;j++) x0[j]=x0[j]/sum;
    if( istep > 2000 )
      {
      printf("Itrative 2000 times! Strike any key to exit.\n");
      getchar();
      exit(1);
      }
    if( err>100.0*EPSILON ) goto aaa;
    free(ik);
    free(x1);
end: return(lamb1) ;
}

void main()
{
int N;
int i,j,rank;
double *eigen;
double *V,*A;
FILE *fp;
FILE *fg;
fp=fopen("input.txt","r");
 fscanf(fp,"%d",&N);
 printf("The Number you put in is: %d\n\n",N);
        A=(double*)calloc(N*N,8);
   for(i=0;i<N;i++)
   {for(j=0;j<N;j++)
       {
        fscanf(fp,"%lf",&A[i*N+j]);
        printf("%12.6f",A[i*N+j]);
       }
   printf("\n");}
eigen=(double *)calloc(N,sizeof(double));
V=(double *)calloc(N*N,sizeof(double));

rank=symeigen( eigen, A, V, N);

printf("Rank=%3d\n",rank);
printf("Eigenvalue:\n");
printf("%14.8lf\n",eigen[N-1]);
printf("Eigenvectors:\n");
for(i=0;i<N;i++)
  {
  printf("%14.8lf",V[i*N+N-1]);
  printf("\n");
  }
getchar();
free(eigen);
free(V);
}