jacobi’s method.cpp

该代码用c++语言实现了Jacobis-Method算法,运行环境为vc6.0

//Jacobi’s Method

//Input
//A symmetric matrix A
//Output
//A diagonal matrix D

/*
Input:
4
0.5 0.25 0 0
0.25 0.8 0.4 0
0 0.4 0.6 0.1
0 0 0.1 1

Output:
0.192242 0 0 0
0 1.18909 0 0
0 0 0.523822 0
0 0 0 0.994844
*/

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

int n;//the size of the matrix A
double A[100][100];//the symmetric matrix A
double P[100][100];//the rotate matrix P
double D[100][100];//the diagonal matrix D
double TOL=0.00001;//the tolerance

//function "sgn"
int sgn(double x)
{
	if(x>0)
		return 1;
	else if(x<0)
		return -1;
	else
		return 0;
}

//P=I
int I()
{
	int j,k;
	
	for(j=1;j<=n;j++)
	{
		for(k=1;k<=n;k++)
			if(k==j)
				P[j][k]=1;
			else
				P[j][k]=0;
	}
	return 1;
}

//P*A*P'
int PAPt()
{
	int i,j,k;
	double PA[100][100];
	double sum;

	//P*A
	for(i=1;i<=n;i++)
	{
		for(j=1;j<=n;j++)
		{
			sum=0;
			for(k=1;k<=n;k++)
				sum=sum+P[i][k]*A[k][j];
			PA[i][j]=sum;
		}
	}
	//P*A*P'
	for(i=1;i<=n;i++)
	{
		for(j=1;j<=n;j++)
		{
			sum=0;
			for(k=1;k<=n;k++)
				sum=sum+PA[i][k]*P[j][k];
			//A=P*A*P'
			A[i][j]=sum;
		}
	}

	return 1;
}

//JacobiMethod
int JacobiMethod()
{
	double b,c;
	double sum;
	int j,k;
	int N;//the max number of iteration
	
	N=0;
	while(N<100)
	{
		//compute the rotate matrix P.
		for(j=2;j<=n;j++)
		{
			for(k=1;k<=j-1;k++)
			{
				I();//P=I
				if(A[j][j]!=A[k][k])
				{
					c=2*A[j][k]*sgn(A[j][j]-A[k][k]);
					b=fabs(A[j][j]-A[k][k]);
					P[j][j]=(sqrt((1+b/sqrt(c*c+b*b))/2));
					P[k][k]=P[j][j];
					//****************跟书上不一样的地方*******************
					//书上是: P[k][j]=c/(2*P[j][j]*sqrt(c*c+b*b));
					//可经过计算和推导，我觉得应该是P[k][j]=-(c/(2*P[j][j]*sqrt(c*c+b*b)));
					P[k][j]=-(c/(2*P[j][j]*sqrt(c*c+b*b)));
					P[j][k]=-P[k][j];
				}
				else
				{
					P[j][j]=sqrt(0.5);
					P[k][k]=P[j][j];
					P[k][j]=P[j][j];
					P[j][k]=-P[k][j];
				}
				
				//compute P*A*P'
				PAPt();
			}
		}

		sum=0;
		for(j=1;j<=n;j++)
		{
			for(k=1;k<=n;k++)
			{
				if(k!=j)
					sum=sum+fabs(A[k][j]);
			}
		}

		//test for success.
		if(sum<TOL)
		{
			cout<<"the diagonal matrix D which approximate the eigenvalues of A."<<endl;
			for(j=1;j<=n;j++)
			{
				for(k=1;k<=n;k++)
				{
					if(k==j)
						D[j][k]=A[j][k];
					else
						D[j][k]=0;
					//output the diagonal matrix D which approximate the eigenvalues of A.
					cout<<D[j][k]<<" ";
				}
				cout<<endl;
			}

			//the method is completed successfully.
			return 1;
		}
		N++;
	}
	
	//the method is completed unsuccessfully.
	return 0;
}

void main()
{
	int i,j;
	
	//Input
	cout<<"please input the size n of the matrix."<<endl;
	cin>>n;
	cout<<"please input a symmetric matrix A."<<endl;
	for(i=1;i<=n;i++)
	{
		for(j=1;j<=n;j++)
			cin>>A[i][j];
	}

	//JacobiMethod
	JacobiMethod();
}