c12.cpp

数值分析最常用的四十种算法

//C12
//Find Eigenvaule,using QR Factorazation

#include<iostream.h>
#include<math.h>

const int N=3;

void MatrixMul(double x[][N],double y[][N],double z[][N],int n)
{
	int i,j,k;
	double sum;
	for(i=0;i<n;i++)
	{
		for (j=0;j<n;j++)
		{
			sum=0;
			for(k=0;k<n;k++)
			{
				sum+=x[i][k]*y[k][j];
			}
            z[i][j]=sum;
		}
	}
}

int sign(double x)
{
	if(x>=0)
		return 1;
	else
		return -1;
}

void QRF(double a[][N],double Q[][N],double R[][N],int n)
{
	double h[N][N],H[N][N],H0[N][N],Q0[N][N],alpha,u[N],po,sum;
	int i,j,k,l;
	for(i=0;i<n;i++)
	{
		for(j=0;j<n;j++)
		{
			if(i==j)
			{
				Q0[i][j]=H0[i][j]=H[i][j]=Q[i][j]=R[i][j]=1;
			}
			else
			{
				Q0[i][j]=H0[i][j]=H[i][j]=Q[i][j]=R[i][j]=0;
			}
		}
	}
	for(i=0;i<n-1;i++)
	{
		for(k=0;k<n;k++)
		{
			for(l=0;l<n;l++)
			{
				if (l==k)
				{
					h[k][l]=1;
				}
				else
				{
					h[k][l]=0;
				}
			}
		}
		sum=0;
		for(j=i;j<n;j++)
		{
			sum+=a[j][i]*a[j][i];
		}
		sum=sqrt(sum);
        alpha=(-1)*sign(a[i][i])*sum;
		for(j=0;j<n-i;j++)
		{
			if(j!=i)
                u[j]=a[j][i];
			else
				u[j]=a[j][i]-alpha;
		}
		po=alpha*(alpha-a[i][i]);
		for(k=i;k<n;k++)
		{
			for(l=i;l<n;l++)
			{
				h[k][l]-=u[k-i]*u[l-i]/po;
			}
		}	
		for(k=0;k<n;k++)
		{
			for(l=0;l<n;l++)
			{
				Q0[k][l]=Q[k][l];
				H0[k][l]=H[k][l];
			}
		}
		//MatrixMul(h,a,H0,n);        Here has an error!!!!!
		MatrixMul(h,H0,H,n);
		//MatrixMul(h,a,H0,n);
		MatrixMul(Q0,h,Q,n);
	}
	MatrixMul(H,a,R,n);
}


void main()
{
	double a[N][N]={5,-3,2,6,-4,4,4,-4,5},
		Q[N][N],R[N][N];
	int i,maxtimes=3;
	for(i=1;i<maxtimes;i++)
	{
		QRF(a,Q,R,N);
		MatrixMul(R,Q,a,N);
	}
	cout<<"The eigenvaule is:"<<endl;
	for(i=0;i<N;i++)
		cout<<a[i][i]<<endl;
}

//失败。原因：第一次所算得矩阵h无法使矩阵a边位列下为零的形式！！