power_cal.c

幂法与反幂法求解矩阵特征值的C语言算法实现:本代码以Hilbert矩阵为计算对象

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define n 11            //矩阵阶数,根据实际修改
#define finished 1
#define stop -1
#define effective_num 9 //有效数字位数
#define MAX_times 100   //最大迭代次数
#define near_eig 0      //近似特征值

typedef struct {
double vec[n];
}_vector;               

typedef struct {
double mat[n][n];
}_matrix;              

static _matrix matrix;
static _matrix L;
static _matrix U;
static _vector vector;
static double threshold;

void create_Hilbert(void)  //生成Hilbert矩阵子函数，根据实际修改
{
	int i,j;
	for(i=0;i<n;i++)
	{
		for(j=0;j<n;j++)
		{
			matrix.mat[i][j]=pow(((i+1)+(j+1)-1),-1);
		}
	}
}

double get_v_max(double x,double y) //按模比较向量元素大小子函数
{
  double max;

  if(fabs(x)>fabs(y))
     max=x;
  else
     max=y;

  y=max;
  return y;
}

void initial(int first_nonzero_index) //初始化，给出初始向量和误差限
{
	int i;
	for(i=0;i<n;i++)
	   vector.vec[i]=0.5;
	threshold=0.5*pow(10, first_nonzero_index-(effective_num));//计算误差限
}

void LU(void) //LU分解子函数
{
    int i,j,k,m;
	double u[n][n];
	double l[n][n];
	double temp_matrix[n][n];
	double temp;

   for(i=0;i<n;i++)
	{
		for(j=0;j<n;j++)
		{
			if(i==j)
			  temp_matrix[i][j]=matrix.mat[i][j]-near_eig;
			else
              temp_matrix[i][j]=matrix.mat[i][j];
		}
	}

    for(i=0;i<n;i++) 
	{
		for(j=0;j<n;j++)
		{
			if(i==j)
                l[i][j]=1.0;    
			if(i<j)
				l[i][j]=0.0;
			if(i>j)
				u[i][j]=0.0;
		}
	}

   for(j=0;j<n;j++)
	   u[0][j]=temp_matrix[0][j];
   for(i=1;i<n;i++)
	   l[i][0]=temp_matrix[i][0]/u[0][0];

   for(k=1;k<n;k++)
   {
	   for(j=k;j<n;j++)
	   {
	     temp=0.0;
	     for(m=0;m<=k-1;m++)
		 {
			 temp=temp+l[k][m]*u[m][j];
		 }
		 u[k][j]=temp_matrix[k][j]-temp;
	   }//生成矩阵U
	   for(i=k+1;i<n;i++)
	   {
		   temp=0.0;
		   for(m=0;m<=k-1;m++)
		   {
			   temp=temp+l[i][m]*u[m][k];
		   }
		   l[i][k]=(temp_matrix[i][k]-temp)/u[k][k];
	   }//生成矩阵L
   }

    for(i=0;i<n;i++) 
	{
		for(j=0;j<n;j++)
		{
	        L.mat[i][j]=l[i][j];
		    U.mat[i][j]=u[i][j];
		}
   	
	}
}

int p_power(void) //幂法子函数
{
	int i,j,k;
    double temp_eig=0.0;
	double max_v,eig;
	double temp_vector[n];
	double Y[n];
	double temp;

	for(k=1;k<=MAX_times;k++)
	{
	   for(i=0;i<n;i++)
	       temp_vector[i]=vector.vec[i];

	   max_v=temp_vector[1];
	   for(i=0;i<n-1;i++)
			max_v=get_v_max(temp_vector[i],max_v);
	   eig=max_v;

		for(i=0;i<n;i++)
			Y[i]=vector.vec[i]/eig;

		for(i=0;i<n;i++)
		{
			temp=0.0;
			for(j=0;j<n;j++)
			{
				temp=temp+matrix.mat[i][j]*Y[j];
			}
			vector.vec[i]=temp;
		}

		if(fabs(eig-temp_eig)<threshold)
		{
			printf("The max eigenvalue of Hilbert Matrix : %.8f.\n",eig);
		    printf("The calculation steps : %d.\n\n",k);
			break;
		}
		else
			temp_eig=eig;
	}
	   
	if(k>MAX_times) //超过允许迭代次数则报错
	{
	   printf("***Faild!!!***\n");
	   return stop;
	   }

   return finished;
}

int n_power(void)  //反幂法子函数
{
	int i,j,k;
	double max_v,eig,alpha,temp;
	double temp_vector[n];
	double temp_eig=1.0;
    _vector X;
    _vector Y;
    _vector Z;

	LU();//作LU三角分解
   
	for(i=0;i<n;i++)
        X.vec[i]=vector.vec[i];
 
	for(k=1;k<=MAX_times;k++)
	{
       for(i=0;i<n;i++)
	       temp_vector[i]=X.vec[i];

	   max_v=temp_vector[1];
	   for(i=0;i<n-1;i++)
			max_v=get_v_max(temp_vector[i],max_v);
	   alpha=max_v;
	
	   for(i=0;i<n;i++)
			Y.vec[i]=X.vec[i]/alpha;	
	  
	   Z.vec[0]=Y.vec[0];
	   for(i=1;i<n;i++)
	   {
		   temp=0.0;
		   for(j=0;j<=i-1;j++)
			   temp=temp+L.mat[i][j]*Z.vec[j];
		   Z.vec[i]=Y.vec[i]-temp;
	   }//解方程组求Z

       X.vec[n-1]=Z.vec[n-1]/U.mat[n-1][n-1];
	   for(i=n-2;i>=0;i--)
	   {
		   temp=0.0;
		   for(j=n-1;j>=i+1;j--)
			   temp=temp+U.mat[i][j]*X.vec[j];
		   X.vec[i]=(Z.vec[i]-temp)/U.mat[i][i];
	   }//解方程组获得X的新值

	   if(fabs(pow(alpha,-1)-pow(temp_eig,-1))<threshold)
	   {
		   eig=near_eig+pow(alpha,-1);
		   printf("The min eigenvalue of Hilbert Matrix : %.23f.\n",eig);
		   printf("The calculation steps : %d.\n\n",k);
		   break;
	   }
	   else
		   temp_eig=alpha;
	}
	
	if(k>MAX_times)  //超过允许迭代次数则报错
	{
	   printf("***Faild!!!***\n");
	   return stop;
	   }

  return finished;
}	

void main(void)
{
  printf("********************************\n");
  printf("*****EIGENVALUE CALCULATION*****\n");
  printf("********************************\n\n");

  printf("Creating Hilbert Matrix...");
  create_Hilbert(); //生成Hilbert矩阵
  printf("Successful!\n\n");

  printf("$CALCULATION RESULT$\n\n");
  initial(1);     
  p_power();
  initial(-14);   
  n_power();     
}