feature.cpp

利用QR分解法求一个10*10矩阵的特征值

#include "stdio.h"
#include "math.h"
#define N 10
double a[N][N],mk[N][N],fv[N][2],x[2][2];
double e=1.0E-12;
double azhen();
double ni();
double qrbijin(int);
double root(double,double);
void main()
{int i, j,k,l,m=N-1;
double g,det,s,t;

 azhen();
 ni();
 while(m>=0)
 {
	 if(m==0)
	 {
		 fv[0][0]=a[0][0];
		 fv[0][1]=0.0;
		 break;
	 }
     else if(m==1)
	 {
		 g=-1.0*(a[m-1][m-1]+a[m][m]);
		 det=a[m-1][m-1]*a[m][m]-a[m][m-1]*a[m-1][m];
		 root(g,det);
		 for(i=0;i<2;i++)
			 for(i=0;i<2;i++)
				 fv[i][j]=x[i][j];
			 break;
	 }
	 else
	 {
		 for(k=1;;k++)
		 {
			 if(fabs(a[m][m-1])<=e)
			 {
				 fv[m][0]=a[m][m];
				 fv[m][1]=0.0;
				 m--;
				 break;
			 }
			 else if(fabs(a[m-1][m-2])<=e)
			 {
				 g=-1.0*(a[m-1][m-1]+a[m][m]);
				 det=a[m-1][m-1]*a[m][m]-a[m][m-1]*a[m-1][m];
				 root(g,det);
				 fv[m][0]=x[1][0];
				 fv[m][1]=x[1][1];
				 fv[m-1][0]=x[0][0];
				 fv[m-1][1]=x[0][1];
				 m-=2;
				 break;
			 }
			 else
			 {
				 s=a[m-1][m-1]+a[m][m];
				 t=a[m-1][m-1]*a[m][m]-a[m][m-1]*a[m-1][m];
				 for(i=0;i<m+1;i++)
					 for(j=0;j<m+1;j++)
						 mk[i][j]=0.0;
				 for(i=0;i<m+1;i++)
					 for(j=0;j<m+1;j++)
					 {
						 for(l=0;l<m+1;l++)
						 {
							 mk[i][j]+=a[i][l]*a[l][j];
						 }
						 mk[i][j]-=s*a[i][j];
						 if(i==j)
							 mk[i][j]+=t;
					 }
				  qrbijin(m+1);
			 }
		 }
	 }
 }					
					                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                            
 for(i=0;i<N;i++)
 {
	 for(j=0;j<N;j++)
     {
		printf("%17.11e\t",a[i][j]);
	}
       printf("\n");
 }
printf("\nFeature Values are:\n");
for(i=0;i<N;i++)
{
	for(j=0;j<2;j++)
		printf("%17.11e\t",fv[i][j]);
	printf("\n");
}
}
/*******************************************************************/
 double azhen()                 /*计算A阵值*/
{int i, j;
    for(i=0;i<N;i++)
    {for(j=0;j<N;j++)
       if(i==j)
         a[i][j]=1.5*cos((i+1)+1.2*(j+1));
	  else
		  a[i][j]=sin(0.5*(i+1)+0.2*(j+1));
    }
	return a[N][N];
  }
/***********************************************************/
double ni()              /*拟上三角化*/
{
int i,j,r,flag;
double d,c,h,t;
double u[N],p[N],q[N],w[N];

for(r=0;r<N-2;r++)
{
	flag=0;
	for(i=r+2;i<N;i++)
		if(fabs(a[i][r])>e)
		{
			flag=1; 
	        break;
		}
		if(flag==0)
			continue;
		d=0.0;
		for(i=r+1;i<N;i++)
			d+=a[i][r]*a[i][r];
		d=sqrt(d);
		if((fabs(a[r+1][r])<=e))   
			c=d;
		else if(a[r+1][r]>e)
			c=-1.0*d;
		else 
			c=d;
		h=c*c-c*a[r+1][r];
		for(i=0;i<N;i++)
			u[i]=p[i]=q[i]=w[i]=0.0;
		u[r+1]=a[r+1][r]-c;
		for(i=r+2;i<N;i++)
			u[i]=a[i][r];
		for(j=0;j<N;j++)
		{
			for(i=0;i<N;i++)
				p[j]+=a[i][j]*u[i];
			p[j]=p[j]/h;
		}
		for(i=0;i<N;i++)
		{
			for(j=0;j<N;j++)
				q[i]+=a[i][j]*u[j];
			q[i]=q[i]/h;
		}
		t=0.0;
		for(i=0;i<N;i++)
			t+=p[i]*u[i];
		t=t/h;
		for(i=0;i<N;i++)
			w[i]=q[i]-t*u[i];
		for(i=0;i<N;i++)
			for(j=0;j<N;j++)
				a[i][j]-=(w[i]*u[j]+u[i]*p[j]);
}
for(i=0;i<N;i++)
	for(j=0;j<N;j++)
		if(fabs(a[i][j])<=e)
			a[i][j]=0.0;
return a[N][N];
}
/*********************************************************/
double qrbijin(int k)                        /*******************对M(k)矩阵QR分解,并求矩阵A(k+1)*******************/
{
	int i,j,r,flag;
    double d,c,h,t;
	double u[N],v[N],p[N],q[N],w[N];
	
	for(r=0;r<k-1;r++)
	{
		flag=0;
        for(i=r+1;i<k;i++)
		{
	      if(fabs(mk[i][r])>e)
		  {
		     flag=1; 
	         break;
		  }
		}
		if(flag==0)
        continue;
		d=0.0;
        for(i=r;i<k;i++)
		{
           d+=mk[i][r]*mk[i][r];
		}
        d=sqrt(d);
        if((fabs(mk[r][r])<=e))   
        c=d;
        else if(mk[r][r]>0)
        c=-1.0*d;
        else 
        c=d;
        h=c*c-c*mk[r][r];
		for(i=0;i<k;i++)
            u[i]=v[i]=p[i]=q[i]=w[i]=0.0;
		u[r]=mk[r][r]-c;
        for(i=r+1;i<k;i++)
            u[i]=mk[i][r];
		for(i=0;i<k;i++)
			for(j=0;j<k;j++)
				v[i]+=mk[j][i]*u[j]/h;
        for(i=0;i<k;i++)
			for(j=0;j<k;j++)
				mk[i][j]-=u[i]*v[j];
        for(j=0;j<k;j++)
		{
			for(i=0;i<k;i++)
				p[j]+=a[i][j]*u[i];
			p[j]=p[j]/h;
		}
		for(i=0;i<k;i++)
		{
			for(j=0;j<k;j++)
				q[i]+=a[i][j]*u[j];
			q[i]=q[i]/h;
		}
		t=0.0;
		for(i=0;i<k;i++)
			t+=p[i]*u[i];
		t=t/h;
		for(i=0;i<k;i++)
			w[i]=q[i]-t*u[i];
		for(i=0;i<k;i++)
			for(j=0;j<k;j++)
				a[i][j]-=(w[i]*u[j]+u[i]*p[j]);
	}
	for(i=0;i<N;i++)
		for(j=0;j<N;j++)
			if(fabs(a[i][j])<=e)
				a[i][j]=0.0;
	return a[N][N];
}
/*******************************************************************/
double root (double b,double c)            /*********求一元二次方程的解************/
{
	double disc;
	disc=b*b-4.0*c;
	if(fabs(disc)<=e)
	{
		x[0][0]=x[1][0]=-0.5*b;
		x[0][1]=x[1][1]=0.0;
	}
	else if(disc>e)
	{
		x[0][0]=-0.5*(b-sqrt(disc));
		x[1][0]=-0.5*(b+sqrt(disc));
		x[0][1]=x[1][1]=0.0;
	}
	else
	{
		x[0][0]=x[1][0]=-0.5*b;
		x[0][1]=0.5*sqrt(-disc);
		x[1][1]=-0.5*sqrt(-disc);
	}
	return x[2][2];
}