jsff.cpp

使用雅戈比方法计算矩阵全部特征值算法的实现

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

#define eeee 0.001
#define MAX  -100000

void inputmatrix(double matrix[][20],int n)
{
	int i,j;
	for(i=0;i<n;i++)
		for(j=0;j<n;j++)
		    scanf("%lf",&matrix[i][j]);
}
bool judge(double matrix[][20],int n)
{
	int i,j;
	for(i=0;i<n;i++)
		for(j=0;j<n;j++)
			if((i!=j)&&fabs(matrix[i][j])>eeee)
				return true;
    return false;
}
double getmaxnum(double matrix[][20],int n,int *p,int *q)
{
	int i,j;
	double max = MAX;
	for(i=0;i<n;i++)
		for(j=0;j<n;j++)
		{
			if(i!=j)
			{
				if(fabs(matrix[i][j])>max)
				{
					max = fabs(matrix[i][j]);
					*p = i;
					*q = j;
				}
			}
		}
    return max;
}
void printre(double matrix[][20],int n)
{
	int i;
	for(i = 0;i<n;i++)
		printf("%lf ",matrix[i][i]);
	cout<<endl;
}
int main()
{
	double matrix[20][20];
	int n,p,q,i;
	double apq,t,c,d,s,t1,t2,temp1,temp2;
	cout<<"please enter the jieshu: ";
	cin>>n;
	cout<<"please enter the matrix"<<endl;
	inputmatrix(matrix,n);
	while(judge(matrix,n))
	{
		apq = getmaxnum(matrix,n,&p,&q);
		s = (matrix[q][q]-matrix[p][p])/(2*apq);
		if(fabs(s)<1e-6)
		{
			t = 1.0;
		}
		else
		{
			t1 =(s+sqrt(s*s+1))*(-1);
			t2 =sqrt(s*s+1)-s;
			if(fabs(t1)<fabs(t2))
			{
				t = t1;
			}
			else t = t2;
		}
		c = 1/sqrt(1+t*t);
		d = t/sqrt(1+t*t);

		for(i=0;i<n;i++)
		{
			if(i!=p&&i!=q)
			{
				temp1 = matrix[p][i];
				temp2 = matrix[q][i];
				matrix[i][p] = matrix[p][i] = c*temp1 - d*temp2;
				matrix[i][q] = matrix[q][i] = c*temp2 + d*temp1;
			}
		}
		matrix[p][p] = matrix[p][p] -t*apq;
		matrix[q][q] = matrix[q][q] +t*apq;
		matrix[p][q] = matrix[q][p] = 0;
	}
	printre(matrix,n);
	return 0;
}