jacobi.cpp

计算矩阵特征值和特征向量的Jacobi算法之程序实现

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <conio.h>
const double EPS=1E-10;  //运算精度

double **D2Array(int row) //动态申请row+1行、row+1列的二维数组
{
	double *Array=new double[(row+1)*(row+1)]; //不使用数组的第0行和第0列
	double **pArray=new double*[row+1];
	for (int i=0;i<=row;i++)
		pArray[i]=Array+i*(row+1);
	return(pArray);
}

void delArray(double **A)  //释放数组占用的内存空间
{
	delete []A[0];
	delete []A;
}
void readFile(double **&A,int &N)
//读入矩阵A及矩阵行数N
{
	FILE *fp;
	double temp;

	fp=fopen(".\\Jacobi.txt","r");
	if (fp==NULL)
	{
		printf("指定位置的文件不存在，请检查！！");
		getch();
		exit(0);
	}

	for (int count=0; ;count++)  //计算矩阵元素个数
	{
		if (feof(fp))
			break;
		fscanf(fp,"%lf",&temp);
	}

	if (count<4)
	{
		printf("矩阵元素不能少于4个，请检查文件！！");
		getch();
		exit(0);
	}

	//计算矩阵行列数，为不超过sqrt(count)的最大整数
	N=(int)sqrt(count);
	A=D2Array(N);  //动态申请矩阵空间

	rewind(fp);
	int r,c;
	for (r=1;r<=N;r++)
		for (c=1;c<=N;c++)
		{
			fscanf(fp,"%lf",&A[r][c]);
			if (feof(fp))
				break;
		}
	fclose(fp);
}

void searchFabsMax(double **A,int N,int &i,int &j)
//在数组A中查找除对角线元素外绝对值最大的元素，返回该元素所在的行号i和列号j
{
	int r,c;
	double max=fabs(A[1][2]);
	i=1; j=2;
	for (r=1;r<=N;r++)
		for (c=1;c<=N;c++)
		{
			if (r==c)  //跳过对角线元素
				continue;
			if (fabs(A[r][c])>max)
			{   max=fabs(A[r][c]);  i=r;  j=c;  }
		}
}

void main( )
{
	double **A;  //实对称矩阵
	int N,r,c;   //N为矩阵大小

	readFile(A,N); //从文件中读取矩阵
	printf("\n初始矩阵为：\n");
	for (r=1;r<=N;r++)
	{
		for (c=1;c<=N;c++)
			printf("%6.2G",A[r][c]);
		printf("\n");
	}

	int count=0; //记录迭代次数
	int flag=1; //设置循环标志，当所有非对角线元素都为0时flag为0
	double **A1=D2Array(N);  //迭代矩阵
	double **V=D2Array(N);   //旋转矩阵
	double **cV=D2Array(N);  //特征向量矩阵
	double **temp=D2Array(N);//临时矩阵
	while (flag==1)
	{
		int i,j;
		searchFabsMax(A,N,i,j); //查找绝对值最大的非对角线元素
//		printf("%6.2g,r=%d,c=%d\n",A[i][j],i,j);  //输出找到的元素

		double fai=0.5*atan(-2*A[i][j]/(A[j][j]-A[i][i]));  //计算φ
//		printf("fai=%lf,sin=%lf,cos=%lf\n\n",fai,sin(fai),cos(fai));
//		getch();

		//计算旋转矩阵
		for (r=1;r<=N;r++)
			for (c=1;c<=N;c++)
				if (r==c) //对角线元素设置为1，否则为0
					V[r][c]=1;
				else
					V[r][c]=0;
		V[i][i]=V[j][j]=cos(fai);
		V[i][j]=-sin(fai);
		V[j][i]=sin(fai);

		//计算特征向量矩阵
		if (count==0) //第1次迭代
		{
			for (r=1;r<=N;r++)
				for (c=1;c<=N;c++)
					cV[r][c]=V[r][c];
		}
		else //第2次及以后的迭代
		{
			for (r=1;r<=N;r++)  //矩阵cV赋值给temp
				for (c=1;c<=N;c++)
					temp[r][c]=cV[r][c];
			for (r=1;r<=N;r++)  //旋转变换矩阵累乘
				for (c=1;c<=N;c++)
				{
					cV[r][c]=0;
					for (int k=1;k<=N;k++)
						cV[r][c]+=temp[r][k]*V[k][c];
				}
		}
		
		//计算新的迭代矩阵A1
		A1[i][i]=A[i][i]*pow(cos(fai),2)+A[j][j]*pow(sin(fai),2)+2*A[i][j]*cos(fai)*sin(fai);
		A1[j][j]=A[i][i]*pow(sin(fai),2)+A[j][j]*pow(cos(fai),2)-2*A[i][j]*cos(fai)*sin(fai);
		for (int s=1;s<=N;s++)  //计算第i行、第i列、第j行及第j列的元素值
		{
			if (s!=i)
				A1[i][s]=A1[s][i]=A[i][s]*cos(fai)+A[j][s]*sin(fai);
			if (s!=j)
				A1[j][s]=A1[s][j]=-A[i][s]*sin(fai)+A[j][s]*cos(fai);
		}
		for (s=1;s<=N;s++) //计算其它元素值
			for (int m=1;m<=N;m++)
				if (m!=i&&m!=j&&s!=i&&s!=j)
					A1[m][s]=A[m][s];
		A1[i][j]=A1[j][i]=0.5*(A[j][j]-A[i][i])*sin(2*fai)+A[i][j]*(pow(cos(fai),2)-pow(sin(fai),2));
		
		for (r=1;r<=N;r++) //处理计算结果，小于EPS则认为等于0
			for (c=1;c<=N;c++)
				if (fabs(A1[r][c])<EPS)
					A1[r][c]=0;

		//更新循环标志
		flag=0;
		for (r=1;r<=N;r++)
			for (c=1;c<=N;c++)
			{
				if (r==c)  //不判断对角线元素
					continue;
				if (fabs(A1[r][c])>=EPS)
					flag=1;
			}
		
		//输出本次迭代结果
//		printf("\n第%d次迭代结果为：\n",++count);
		for (r=1;r<=N;r++)
		{
			for (c=1;c<=N;c++)
			{
				A[r][c]=A1[r][c]; //同时把A1赋值给A
				//printf("%15.6lG",A[r][c]);
			}
			//printf("\n");
		}
	} //while循环结束

	//输出特征值
	getch();
	printf("\n特征值为：\n");
	for (r=1;r<=N;r++)
		printf("%10lG",A[r][r]);
	printf("\n");

	//输出特征向量
	printf("\n特征向量为：\n");
	for (r=1;r<=N;r++)
	{
		for (c=1;c<=N;c++)
		{
			if (fabs(cV[r][c])<EPS) //处理计算结果
				cV[r][c]=0;
			printf("%10lG",cV[r][c]);
		}
		printf("\n");
	}
	printf("\n");

	//释放数组占用的内存空间
	delArray(A);
	delArray(A1);
	delArray(V);
	delArray(cV);
	delArray(temp);
}