finalpro.cpp

用多种算法实现任给实对称矩阵的所有特征值

#include<iostream.h>
#include<stdio.h>
#include<stdlib.h>
#include<time.h>
#include<math.h>
void fun(double a[],int n);//排序函数将求得的特征值从小到大排列
int jaccobi(double a[],double u[],int n0,double eps);//普通jaccobi法
void jacobbigg(double a[],double u[],int n0,double eps);//jaccobi过关法
void hhtr(double a[],int n);//矩阵变换函数：用householder相似变化将一般是矩阵化为上Hessenberg阵
int qr(double h[],int n0,double u[],double v[],double eps,int jt);//双步QR法求特征值
void fun(double a[],int n)
{	
	int i,j,d;double m,s;
	for(j=0;j<n;j++)
	{
		s=a[j];
		for(i=j;i<n;i++)
		{			
			if(a[i]<=s)
			{s=a[i];d=i;}
		}
		m=a[j];
		a[j]=a[d];
		a[d]=m;
		}
	for(i=0;i<n;i++)
		cout<<a[i]<<endl;
	cout<<endl;
}
int jacobi(double a[],double u[],int n0,double eps)
{
	int i,j,p,q,pj,qj,ip,iq;
	double m,n,w,st,s,c,temp;
	double t;//存储非对角线的最大值Apq
	p=0;q=0;
	while(1)
	{
		t=0.0;
		for(i=0;i<n0;i++)//选出非对角线元素中绝对值最大者Apq
		{
			for(j=i+1;j<n0;j++)
			{
				if(fabs(a[i*n0+j])>t){t=fabs(a[i*n0+j]);p=i;q=j;}
			}
		}
		if(fabs(a[p*n0+q])<eps)//满足精度则打印结果
		{			
			for(i=0;i<n0;i++)u[i]=a[i*n0+i];
			cout<<"----------------"<<endl;
			fun(u,n0);
			return 1;
		}
		else//否则按照jacobi算法计算矩阵A的新元素
		{
			m=-a[p*n0+q];
			n=0.5*(a[q*n0+q]-a[p*n0+p]);
			if(n>=0)w=m/sqrt(m*m+n*n);
			else w=-m/sqrt(m*m+n*n);
			st=w;
			s=w/sqrt(2*(1+sqrt(1-w*w)));
			c=sqrt(1-s*s);
			temp=a[p*n0+p];
			a[p*n0+p]=temp*c*c+a[q*n0+q]*s*s+a[p*n0+q]*st;
			a[q*n0+q]=temp*s*s+a[q*n0+q]*c*c-a[p*n0+q]*st;
			a[p*n0+q]=0.0;a[q*n0+p]=0.0;
			for(j=0;j<n0;j++)
			{
				if((j!=p)&&(j!=q))
				{
					pj=p*n0+j;qj=q*n0+j;temp=a[pj];
					a[pj]=temp*c+a[qj]*s;
					a[qj]=-temp*s+a[qj]*c;
				}
			}
			for(i=0;i<n0;i++)
			{
				if((i!=p)&&(i!=q))
				{
					ip=i*n0+p;iq=i*n0+q;temp=a[ip];
					a[ip]=temp*c+a[iq]*s;
					a[iq]=-temp*s+a[iq]*c;
				}
			}
	}
}
	return 1;
}
void jacobigg(double a[],double u[],int n0,double eps)
{
	int i,j,p,q,pj,qj,ip,iq;
	double m,n,w,st,s,c,temp;
	double t,sum,r;//存储最大值Apq
	p=0;q=0;
	sum=0.0;
	t=0.0;
	for(i=0;i<n0;i++)
	{
		for(j=i+1;j<n0;j++)
		{
			sum=sum+a[i*n0+j]*a[i*n0+j];
		}
	}
	sum=sqrt(2*sum);r=sum;
	r=r/n0;
	while(r>=eps)
	{
loop1:
	for(i=0;i<n0;i++)
	for(j=i+1;j<n0;j++)
	{
		if(fabs(a[i*n0+j])>=r)
		{
			p=i;q=j;
			m=-a[p*n0+q];
			n=0.5*(a[q*n0+q]-a[p*n0+p]);
			if(n>=0)w=m/sqrt(m*m+n*n);
			else w=-m/sqrt(m*m+n*n);
			st=w;
			s=w/sqrt(2*(1+sqrt(1-w*w)));
			c=sqrt(1-s*s);
			temp=a[p*n0+p];
			a[p*n0+p]=temp*c*c+a[q*n0+q]*s*s+a[p*n0+q]*st;
			a[q*n0+q]=temp*s*s+a[q*n0+q]*c*c-a[p*n0+q]*st;
			a[p*n0+q]=0.0;a[q*n0+p]=0.0;
			for(j=0;j<n0;j++)
			{
				if((j!=p)&&(j!=q))
				{
					pj=p*n0+j;qj=q*n0+j;temp=a[pj];
					a[pj]=temp*c+a[qj]*s;
					a[qj]=-temp*s+a[qj]*c;
				}
			}
			for(i=0;i<n0;i++)
			{
				if((i!=p)&&(i!=q))
				{
					ip=i*n0+p;iq=i*n0+q;temp=a[ip];
					a[ip]=temp*c+a[iq]*s;
					a[iq]=-temp*s+a[iq]*c;
				}
			}
			goto loop1;
		}
	}
	r=r/n0;
	}
	for(i=0;i<n0;i++)u[i]=a[i*n0+i];
	fun(u,n0);
	return;
}
void hhtr(double a[],int n)
{
	int i,j,k,m;
	double f,g,s,tol,t;
	tol=1.0e-12;
	for(k=1;k<=n-2;k++)
	{
		s=0.0;
		for(m=k+1;m<=n;m++)
		{
			s=s+a[(m-1)*n+k-1]*a[(m-1)*n+k-1];
		}
		if(s<=tol)
			g=0.0;
		else
		{
			if(a[k*n+k-1]>=0)g=sqrt(s);
			if(a[k*n+k-1]<0)g=-sqrt(s);
			t=s+g*a[k*n+k-1];
			a[k*n+k-1]=a[k*n+k-1]+g;
			for(j=k+1;j<=n;j++)
			{
				f=0.0;
				for(m=k+1;m<=n;m++)
				{
					f=f+a[(m-1)*n+k-1]*a[(m-1)*n+j-1];
				}
				for(i=k+1;i<=n;i++)
				{
					a[(i-1)*n+j-1]=a[(i-1)*n+j-1]-f*a[(i-1)*n+k-1]/t;
				}
			}
			for(i=1;i<=n;i++)
			{
				f=0.0;
				for(m=k+1;m<=n;m++)f=f+a[(m-1)*n+k-1]*a[(i-1)*n+m-1];
				for(j=k+1;j<=n;j++)a[(i-1)*n+j-1]=a[(i-1)*n+j-1]-f*a[(j-1)*n+k-1]/t;
			}
		}
		a[k*n+k-1]=-g;
	}
	for(m=3;m<=n;m++)
	{
		for(k=1;k<=m-2;k++)
			a[(m-1)*n+k-1]=0;
	}
}
int qr(double h[],int n0,double u[],double v[],double eps,int jt)
{
	int n,kfail,its;//收缩后的矩阵阶数,失败标记,迭代次数
	double d,x,y,w,b,c,z,s,p,q,r,e,f,a;
	n=n0;kfail=0;
	int i,j,k,l,i1,j1;
	its=0;
	while(n!=0)
	{	
		l=n-1;
		while((fabs(h[l*n0+l-1])>(eps*
			(fabs(h[(l-1)*n0+l-1])+fabs(h[l*n0+l]))))&&(l>0))l=l--;
		if(l==n-1)
		{u[n-1]=h[(n-1)*n0+n-1];v[n-1]=0.0;n--;its=0;}
		else
		{
			if(l==n-2)
			{
				x=h[(n-1)*n0+n-1];y=h[(n-2)*n0+n-2];w=h[(n-1)*n0+n-2]*h[(n-2)*n0+n-1];
				a=0.5*(y-x);b=a*a+w;c=sqrt(fabs(b));x=x;
				if(b>0){if(a>=0)z=a+c;if(a<0)z=a-c;u[n-2]=x+z;v[n-2]=0.0;u[n-1]=x-w/z;v[n-1]=0.0;}
				if(b<=0){u[n-1]=x+a;v[n-1]=c;u[n-2]=x+a;v[n-2]=-c;}
				n=n-2;its=0;
			}
			else
			{
				if(its==jt){cout<<"fail"<<endl;kfail=-1;return -1;}
				else
				{
					its++;
					for(j=l+2;j<=n-1;j++)h[j*n0+j-2]=0.0;
					for(j=l+3;j<=n-1;j++)h[j*n0+j-3]=0.0;
					for(k=l;k<=n-2;k++)
					{
						if(k!=l)
						{
							p=h[k*n0+k-1];q=h[(k+1)*n0+k-1];r=0.0;
							if(k!=n-2)r=h[(k+2)*n0+k-1];
						}
						else
						{
							x=h[(n-1)*n0+n-1]+h[(n-2)*n0+n-2];
							y=h[(n-1)*n0+n-1]*h[(n-2)*n0+n-2]-h[(n-1)*n0+n-2]*h[(n-2)*n0+n-1];
							p=h[l*n0+l]*(h[l*n0+l]-x)+h[l*n0+l+1]*h[(l+1)*n0+l]+y;
							q=h[(l+1)*n0+l]*(h[l*n0+l]+h[(l+1)*n0+l+1]-x);
							r=h[(l+1)*n0+l]*h[(l+2)*n0+l+1];
						}
						if((fabs(p)+fabs(q)+fabs(r))!=0.0)
						{
							if(p>=0.0)s=sqrt(p*p+q*q+r*r);
							if(p<0.0)s=-sqrt(p*p+q*q+r*r);
							if(k!=l)h[k*n0+k-1]=-s;
							d=-q/s; e=-r/s; x=-p/s;
							y=-x-e*r/(p+s);
							f=d*r/(p+s);
							z=-x-d*q/(p+s);
							
							for(j=k;j<=n-1;j++)
							{
								i1=k*n0+j; j1=(k+1)*n0+j;
								p=x*h[i1]+d*h[j1];
								q=d*h[i1]+y*h[j1];
								r=e*h[i1]+f*h[j1];
								if(k!=n-2)
								{
									p=p+e*h[(k+2)*n0+j];
									q=q+f*h[(k+2)*n0+j];
									h[(k+2)*n0+j]=r+z*h[(k+2)*n0+j];
								}
								h[(k+1)*n0+j]=q; h[k*n0+j]=p;
							}
							j=k+3;
							if(j>=n-1)j=n-1;
							for(i=l;i<=j;i++)
							{
								i1=i*n0+k; j1=i*n0+k+1;
								p=x*h[i1]+d*h[j1];
								q=d*h[i1]+y*h[j1];
								r=e*h[i1]+f*h[j1];
								
								if(k!=n-2)
								{
									
									p=p+e*h[i*n0+k+2];
									q=q+f*h[i*n0+k+2];
									h[i*n0+k+2]=r+z*h[i*n0+k+2];
								}	
								h[i*n0+k+1]=q; h[i*n0+k]=p;
							}
						}
					}
				}
			}
		}
	}
	return 1;
}
int main()
{
	int i,j,n,jt;
	clock_t ts,te;
	char c,*filename,*method;
	static double u[200],v[200];
	double a[40000];
	double eps;
	double time;
	FILE *fp;
	cout<<"----------------任给实对称矩阵求其所有特征值------------------"<<endl;
	cout<<"按任意键继续..."<<endl;
	getchar();
	system("cls");
loop:
	cout<<"请输入所求矩阵的结束：1. 30阶   2. 50阶   3. 100阶   4. 200阶   5. 退出"<<endl;
	cin>>i;
	switch(i)
	{
	case 1:filename="matrix_30.txt";break;
	case 2:filename="matrix_50.txt";break;
	case 3:filename="matrix_100.txt";break;
	case 4:filename="matrix_200.txt";break;
	case 5:exit(-1);
	default:break;
	}
	if((fp=fopen(filename,"r+"))==0)
	{
		cout<<"cannot open the file!"<<endl;//报错
		getchar();
		return 0;
	}
	fscanf(fp,"%d",&n);
	for(i=0;i<n;i++)//文件读入，用双精度型读入
	{
		for(j=0;j<n;j++)
			fscanf(fp,"%lf",a+i*n+j);
		c=NULL;
		while((c!=';')&&(!feof(fp)))//从一行内最后一个数据一直读取字符，知道遇到“；”或者碰到文件尾部结束，进行下一行的读取或者结束。
			c=fgetc(fp);
	}
	fclose(fp);
	cout<<endl;
	system("cls");
	cout<<"请输入想采用的求特征值方法："<<endl;
	i=0;
	while(i!=4)
	{
	cout<<"1.jacobi过关法   2.jacobi法     3.双步QR法   4.返回上一级菜单"<<endl;
	cin>>i;
	switch(i)
	{
	case 1:
		cout<<"希望的精度是："<<endl;cin>>eps;
		ts=clock();
		jacobigg(a,u,n,eps);
		te=clock();
		time=(double)(te-ts)/CLOCKS_PER_SEC;
		cout<<"耗时为："<<time<<"s"<<endl;
		method="雅克比过关法";
		break;
	case 2:
		cout<<"希望的精度是："<<endl;cin>>eps;
		ts=clock();
		j=jacobi(a,u,n,eps);
		te=clock();
		time=(double)(te-ts)/CLOCKS_PER_SEC;
		cout<<"耗时为："<<time<<"s"<<endl;
		method="普通雅克比法";
		break;
	case 3:
		cout<<"希望的迭代次数不超过："<<endl;cin>>jt;
		cout<<"希望的精度是："<<endl;cin>>eps;
		ts=clock();
		hhtr(a,n);
		j=qr(a,n,u,v,eps,jt);
		te=clock();
		time=(double)(te-ts)/CLOCKS_PER_SEC;
		if (j>0)
		fun(u,n);
		cout<<"耗时为:"<<time<<"s"<<endl;
		method="双步QR法";
		break;
	case 4:
		system("cls");
		goto loop;break;
	default: break;
	}
	cout<<"按任意键继续..."<<endl;
	getchar();
	system("cls");
	if(n==30)filename="char_30.txt";
	if(n==50)filename="char_50.txt";
	if(n==100)filename="char_100.txt";
	if(n==200)filename="char_200.txt";		
	if((fp=fopen(filename,"w+"))==NULL)
	{cout<<"cannot open this file!"<<endl;exit(0);}
	for(i=0;i<n;i++)
		fprintf(fp,"%lf\n",u[i]);
	fprintf(fp,"%s","矩阵阶数");
	fprintf(fp,"%d\n",n);
	fprintf(fp,"%s\n",method);
	fprintf(fp,"%s\n","运行时间为:");
	fprintf(fp,"%lf",time);
	fprintf(fp,"%s\n","s");
	}
	
		return 0;
	}