main.c

GivensHouseholder法求解特征向量和特征值

/*bisect.c:Givens-Householder mothed for eigen-problems of symmetric*/
#include <stdio.h>
#include <stdlib.h>
#include <conio.h>
#include <time.h>
#include <math.h>

#define EPSILON 1.0e-10

#define N 5
#define TM 20

void avector( double A[N][N], double V[N][N] );
void givens( double A[N][N] );
void bisection( double ai[N], double bi[N], double cent[N] );
int samesign( double lembda, double ai[N], double bi2[N] );
void tvector( double lembda, double ai[N], double bi[N], double x[N]);
void bisect( double A[N][N], double ai[N], double bi[N], double cent[N] )
{
int i,j;
givens( A );
for(i=0;i<N;i++)
  {
  for(j=0;j<N;j++) printf("%12.8f",A[i][j]);
  printf("\n");
  }

for(i=0;i<N-1;i++)
  {
  ai[i]=A[i][i];
  bi[i]=A[i][i+1];
  }
ai[N-1]=A[N-1][N-1];
bi[N-1]=0.0;
bisection( ai, bi, cent );
}
/********************************************************************/
void givens( double A[N][N] )
{
int i,l,k;
double d,c,s,temp;
for(k=0;k<N-2;k++)
  {
  for(l=k+2;l<N;l++)
    {
    if( A[l][k] == 0.0 ) A[k][l]=1.0;
    else
      {
      d=A[k+1][k]*A[k+1][k]+A[l][k]*A[l][k];
      d=sqrt(d);
      c=A[k+1][k]/d;
      s=A[l][k]/d;
      A[k][l]=c;
      A[l][k]=s;
      A[k+1][k]=d; A[k][k+1]=d;
      for(i=k+1;i<N;i++)
        {
        temp=A[i][k+1];
        A[i][k+1]= temp*c+A[i][l]*s;
        A[i][l] = -temp*s+A[i][l]*c;
        }
      for(i=k+1;i<N;i++)
        {
        temp=A[k+1][i];
        A[k+1][i]= temp*c+A[l][i]*s;
        A[l][i] = -temp*s+A[l][i]*c;
        }
      }    /* end of if */
    }      /* end of l  */
  }        /* end of k  */
}
/*******************************************************************/
void bisection( double ai[N], double bi[N], double cent[N] )
{
int i,j,l,k,t,sl[N];
static double bi2[N],ll[N],ur[N],norm,temp,c;
for(i=0;i<N-1;i++)
  {
  if(fabs(bi[i]) < EPSILON )
     {
     printf("The b[%d] is too small absolutely!\n",i);
     printf("Strike any key to exit.\n");
      exit(1);
     }
  }
bi2[0]=bi[0]*bi[0];
norm=fabs(ai[0])+fabs(bi[0]);
for(i=1;i<N;i++)
  {
  bi2[i]=bi[i]*bi[i];
  temp=fabs(bi[i-1])+fabs(ai[i])+fabs(bi[i]);
  if( norm < temp ) norm=temp;
  }
for(i=0;i<N;i++)
  {
  ll[i]=-norm;
  ur[i]=norm;
  }
c=ll[0];
sl[0]=samesign( c, ai, bi2 );
for(t=0;t<TM;t++)
  {
  for(i=0;i<N;i++)
    {
    aa: c=0.5*(ll[i]+ur[i]);
    k=samesign( c, ai, bi2 );
    if( k == sl[i] ) { ll[i]=c; goto aa; }
    else ur[i]=c;
    if( k < (sl[i]-1) ) goto aa;
    if( (i < N-1) && ( ll[i+1] < ur[i]) ) ll[i+1]=ur[i];
    sl[i+1]=k;
    }       /* end for  i loop */
  }         /* end for  t loop */
for(i=0;i<N;i++) cent[i]=0.5*(ll[i]+ur[i]);
}
/********************************************************************/
int samesign( double lembda, double ai[N], double bi2[N] )
{
int i, slembda;
static double si[N+1],si2;
si[0]=1.0;
si[1]=ai[0]-lembda;
for(i=2;i<=N;i++)
  {
  si2=fabs(si[i-1]*si[i-2]);
  if( si2 > EPSILON*EPSILON ) si[i]=ai[i-1]-lembda-bi2[i-2]/si[i-1];
  if( fabs(si[i-2]) < EPSILON ) si[i]=ai[i-1]-lembda;
  if( fabs(si[i-1]) < EPSILON ) si[i]=-EPSILON;
  }
slembda=0;
for(i=1;i<=N;i++)
  if( si[i] >= 0.0 ) slembda=slembda+1;
return(slembda);
}
/*******************************************************************/
void tvector( double lembda, double ai[N], double bi[N], double x[N])
{
int i;
double norm;
for(i=0;i<N-1;i++)
  {
  if(fabs(bi[i]) < EPSILON )
     {
     printf("The b[%d] is too small absolutely!\n",i);
     printf("Strike any key to exit.\n");
      exit(1);
     }
  }
x[0]=1.0;
x[1]=-(ai[0]-lembda)*x[0]/bi[0];
for(i=1;i<N-1;i++) x[i+1]=-(bi[i-1]*x[i-1]+(ai[i]-lembda)*x[i])/bi[i];
norm=0.0;
for(i=0;i<N;i++) norm=norm+x[i]*x[i];
norm=sqrt(norm);
for(i=0;i<N;i++) x[i]=x[i]/norm;
}
/********************************************************************/
void avector( double A[N][N], double V[N][N] )
{
int i,l,k;
double c,s,temp;
for(k=N-3;k>=0;k--)
  {
  for(l=N-1;l>=k+2;l--)
    {
    c=A[k][l];
    s=A[l][k];
    for(i=0;i<N;i++)
      {
      temp=V[k+1][i];
      V[k+1][i]=temp*c-V[l][i]*s;
      V[l][i] = temp*s+V[l][i]*c;
      }
    }      /* end of l  */
  }        /* end of k  */
}

/*givenshm.c:Givens-Householder mothed for eigen-problems of symmetric*/

void main( void )
{
	int i,j,k;
	static double A[N][N] = { { 10.0,  1.0,  2.0,  3.0,  4.0 },
	{  1.0,  9.0, -1.0,  2.0, -3.0 },
	{  2.0, -1.0,  7.0,  3.0, -5.0 },
	{  3.0,  2.0,  3.0, 12.0, -1.0 },
	{  4.0, -3.0, -5.0, -1.0, 15.0 } };

	static double ai[N],bi[N],V[N][N],eigen[N],x[N],lembda;

	printf("Givens-Householder method for eigenvalue problems:\n");
	printf("Givens transform:\n");
	bisect( A, ai, bi, eigen );
	printf("Eigenvalues by using bisection:\n");
	for(i=0;i<N;i++) printf("%4d %15.9f\n",i,eigen[i]);

	for(k=0;k<N;k++)
	{
		lembda=eigen[k];
		tvector( lembda, ai, bi, x );
		for(i=0;i<N;i++) V[i][k]=x[i];
	}
	printf("Eigenvectors of matrix T:\n");
	for(i=0;i<N;i++)
	{
		for(j=0;j<N;j++) printf("%12.8f",V[i][j]);
		printf("\n");
	}
	avector( A, V );
	printf("Eigenvectors of matrix A:\n");
	for(i=0;i<N;i++)
	{
		for(j=0;j<N;j++) printf("%12.8f",V[i][j]);
		printf("\n");
	}
	
}