eigenvalue.cpp

vc＋＋实现神经网络的LMS算法

/*幂法求实对称矩阵特征值*/
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <time.h>
#include <string.h>
#include <malloc.h>
#include "conio.h"
#include "eigen.h"

#define EPSILON 1.0e-10

float lamb1;

/*****************************************************************
入口参数：eigen（N个特征值），A(输入的实对称矩阵)
返回参数：矩阵的秩
******************************************************************/
int symeigen(float *eigen, float **A, int N)
{
	int i,j,step,rank;
	static float **B =(float **)malloc(4*N*N),
	static float *x0 =(float *)malloc(4*N);
	static float mlamb;

	for(i=0; i<N; i++)
	 { 
		*( eigen+i )=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++) 
			{
				 *(float *)(B+i*N+j) = *(float *)(A+i*N+j);
			}

			*(float *)(B+i*N+i) = *(float *)(A+i*N+i)-lamb1;
		}

		  lamb1 = invpower( 1, x0, (float **)B, N);

		  if( fabs(lamb1) > EPSILON ) 
		  {
			  *( eigen+step )= mlamb+ (1/lamb1);
		  }

		  if(lamb1==0.0) 
		  {
			  *( eigen+step )=mlamb;
		  }

		bbb:
		  if( fabs( *( eigen+step ) ) > EPSILON ) 
		  {
			  rank=rank+1;
		  }
		  
	}

	return rank;
}
/************************************************************************/
void power(int i0, float *x0, float **A, int N)
{
	static float lamb,sum,err,err1,eps;
	static float *x1 =(float *)malloc(4*N);
	float *V =(float *)malloc(4*N*N);
	int i,j,jj,ii;

	memset(V, 0, 4*N*N);
	
	eps=(float)EPSILON*1000;
	
	raa:
		ii=0;
		srand(clock());                                         /*for creat x0[]*/
	
		for(j=0; j<i0; j++) 
		{ 
			for(i=0; i<N; i++) 
			{
				*x0=(float)rand();
			}
		}

		sum=0.0;
		for(i=0; i<N; i++)
		{ 
			*( x0+i )=(float)rand(); 
			sum +=( ( *(x0+i) )* ( *(x0+i) ) ); 
		}

		sum=(float)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 ) += (  ( *(float *)( A+j*N+jj ) )* ( *( x0 +jj ) )  );
			  }
				  
		  }
		lamb1=0.0;
		for(j=0; j<N; j++)
		{
			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=(float)fabs(lamb1-lamb);
			err1=(float)err/(float)fabs(lamb1); 
		}
		
		lamb= lamb1;

	bbb:
		if(i0>0)
		{
			for(i=0; i<i0; i++)
			{
				sum=0.0;
				for(j=0;j<N;j++) 
				{
					sum += ( (*(x1+j)) * (*(float *)( V+j*N+i )) );
				}

				for(j=0;j<N;j++) 
				{
					*(x1+j) -= ( sum* (*(float *)( V+j*N+i )) );
				}
			}
		  }

	sum=0.0;
	for(j=0; j<N; j++) 
	{
		sum +=( (*(x1+j)) * (*(x1+j)) );
	}

	sum=(float)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;
	}
	
	endd:;
}


/***************************************************************/
float invpower( int init, float *x0, float **B, int N)
{
	int i,j,k,i0,istep;
	int *ik =(int *)malloc(4*N);
	static float lamb,lamb1,sum,err,temp;
	static float  *x1 =(float *)malloc(4*N);

	memset(ik, 0, 4*N);
	
	istep=0;

	if( init == 0 )              // for creat x0[]
	{
		srand(clock());
		raa: 
		  sum=0.0;
		  for(i=0; i<N; i++)
		  { 
			  *(x0+i) = (float)rand(); 
			  sum +=( ( *(x0+i) ) * ( *(x0+i) ) );
		  }

		  sum=(float)sqrt(sum);
		  if( sum < EPSILON ) 
		  {
			  goto raa;
		  }

		  for(i=0; i<N; i++) 
		  {
			  *(x0+i) /=sum;      // x0[] is created
		  }
	}
	for(k=0; k<N; k++)                 //colum pivot Doolittle decompsition
	{
		temp=(float)fabs( *(float *)(B+k*N+k) );       
		*(ik+k) = k;
		for(i=k; i<N; i++)
		{
			if( fabs( *(float *)(B+i*N+k)) > temp ) 
			{ 
				temp=(float)fabs( *(float *)(B+i*N+k) );
				*(ik+k) = i; 
			}
		  i0 = *(ik+k);
		  if( i0 != k )
		  {
				for(j=0;j<N;j++)
				 { 
					temp= *(float *)(B+k*N+j); 
					*(float *)(B+k*N+j)= *(float *)(B+i0*N+j); 
					*(float *)(B+i0*N+j) =temp; 
				}
		  }

		  if( fabs(*(float *)(B+k*N+k)) < EPSILON ) 
		  { 
			  lamb1=0.0; 
			  goto end; 
		  }

		  for(i=k+1; i<N; i++)
		  {
				*(float *)(B+i*N+k) /=*(float *)(B+k*N+k);
				for(j=k+1;j<N;j++) 
				{
					*(float *)(B+i*N+j) -= ( ( *(float *)(B+i*N+k) ) * ( *(float *)(B+k*N+j) ) );
				}
		  }
		}			  // decompsition is ended
	aaa:
		istep ++;
		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] -= ( ( *(float *)(B+i*N+j) ) * ( *(x0+j) ) );
		  }
		}
		for(i=N-1; i>=0; i--)
		{
			  for(j=i+1; j<N; j++) 
			  {
				  x0[i] -=( ( *(float *)(B+i*N+j) ) * ( *(x0+j) ) );
			  }
			  x0[i] /= *(float *)(B+i*N+i);
		}                                         // end for solution 
		
		lamb1=0.0;
		for(j=0; j<N; j++) 
		{
			lamb1 += ( ( *(x0+j) ) * ( *(x1+j) ) );
		}

		err=(float)fabs(lamb1-lamb)/((float)fabs(lamb1)+1);
		sum=0.0;

		for(j=0; j<N; j++) 
		{
			sum += ( ( *(x0+j) )*( *(x0+j) ) );
		}

		sum=(float)sqrt(sum);
		
		for(j=0; j<N; j++) 
		{
			*(x0+j) /=sum;
		}

		if( istep > 2000 )
		{
		  printf("Itrative 2000 times! Strike any key to exit.\n");
		  getch();
		  exit(1);
		}

		if( err>1000.0*EPSILON ) 
		{
			goto aaa;
		}
	}

	end: 
		return lamb1 ;
}