matrixl.cpp

此程序为求特征值与特征向量的源码程序

// matrixl.cpp : Defines the entry point for the console application.
//

#include "stdafx.h"
#include "math.h"


int main(int argc, char* argv[])
{
	int Inv(double **a, int n, double **b, double **c);
	int Eigs(double **a, int n, double *d, double **v, int& nrot);
//	int Eigs(double **a, int n, double **v, double eps, int jt);
	int i=0;
	int j=0;
	int k=0;
	double **a=new double*[3], **b=new double*[3], **c=new double*[3];
	double **s=new double*[3], **vector=new double*[3];
	double **a1=new double*[3];
	for (i=0; i<3; i++)
	{
		a[i]=new double[3];
		b[i]=new double[3];
		c[i]=new double[3];
		a1[i]=new double[3];
		s[i]=new double[3];
		vector[i]=new double[3];
		for (j=0; j<3; j++)
		{
			a[i][j]=0.0;
			b[i][j]=0.0;
			c[i][j]=0.0;
			a1[i][j]=0.0;
			s[i][j]=0.0;
			vector[i][j]=0.0;
		}
	}
	a[0][0]=2;
	a[1][1]=1;
	a[2][2]=2;
	b[0][0]=1;
	b[1][1]=1;
	b[2][2]=1;
	c[0][0]=1.0;
	c[0][1]=2.0;
	c[0][2]=3.0;
	c[1][0]=2.0;
	c[1][1]=2.0;
	c[1][2]=3.0;
	c[2][0]=3.0;
	c[2][1]=3.0;
	c[2][2]=3.0;
	for (i=0; i<3; i++)
		for (j=0; j<3; j++)
			for (k=0; k<3; k++)
				s[i][j]+=a[i][k]*c[k][j];
	int value=0;
	double *d=new double[3];
	int nrot;
	value=Inv(a, 3, b, a1);
	Eigs(c, 3, d, vector, nrot);
//	Eigs(c, 3, vector, 0.1, 30);
	if (value==1)
	{
	for (i=0; i<3; i++)
		for (j=0; j<3; j++)
			printf("a[%d][%d]=%f\tc[%d][%d]=%f\ts[%d][%d]=%f\na1[%d][%d]=%f\tvector[%d][%d]=%f\n", 
			i,j,a[i][j], i,j,c[i][j], i,j,s[i][j], i,j,a1[i][j], i,j,vector[i][j]);
	}
	printf("特征向量为：%f, %f, %f\n", d[0], d[1], d[2]);
	return 1;
}

/********************************************************************************
* 求矩阵的逆
* 返回值小于0表示矩阵不可逆
* 返回值大于0表示正常返回
* a-长度为n*n的数组，存放实矩阵
* n-矩阵的阶数
* b-长度为n*n的数组，存放单位矩阵
* c-长度为n*n的数组，存放实矩阵的逆矩阵 
*********************************************************************************/
int Inv(double **a, int n, double **b, double **c)
{
	int i,j,m,l,k,q;
	double x=0.0,y=0.0;
	for(i=0;i<n;i++)
		for(j=0;j<n;j++)
		{
			if(i==j)
				b[i][j]=1;
			else
				b[i][j]=0;
		}
	for(i=0;i<n-1;i++)
	{
		if(a[i][i]==0)
			for(l=1;l<n-i;l++)
				if(a[i+l][i]!=0)
					for(m=0;m<n;m++)
					{
						x=a[i+l][m];a[i+l][m]=a[i][m];a[i][m]=x;
						y=b[i+l][m];b[i+l][m]=b[i][m];b[i][m]=y;
						break;
					}
					for(j=i+1;j<n;j++)
					{
						for(k=i;k<n;k++)
							a[j][k]=a[j][k]-a[i][k]*a[j][i]/a[i][i]; 
						for(q=0;q<n;q++)
							b[j][q]=b[j][q]-a[j][i]*b[i][q]/a[i][i];
					}
	}
	if(a[n-1][n-1]==0)
		return(-1);
	else
		for(i=n-1;i>0;i--)
			for(j=i-1;j>=0;j--)
				for(k=0;k<n;k++)
					b[j][k]=b[j][k]-b[i][k]*a[j][i]/a[i][i] ;

	for(i=0;i<n;i++)
		for(j=0;j<n;j++)
			c[i][j]=b[i][j]/a[i][i];
	return(1);
}

/********************************************************************************
* 求实对称矩阵的特征值及特征向量的雅格比法
* 利用雅格比(Jacobi)方法求实对称矩阵的全部特征值及特征向量
* 返回值小于0表示超过迭代jt次仍未达到精度要求
* 返回值大于0表示正常返回
* a-长度为n*n的数组，存放实对称矩阵，返回时对角线存放n个特征值
* n-矩阵的阶数
* v-长度为n*n的数组，返回特征向量(按列存储)
* eps-控制精度要求
* jt-整型变量，控制最大迭代次数
*********************************************************************************/
int Eigs(double **a, int n, double *d, double **v, int& nrot)
{
	double b[100], z[100], g, s, c, h, tau, sss, ddd, t, tresh, sm, theta;
	int ip, iq, i, j;
	for (ip=0; ip<n; ip++)
	{
		for (iq=0; iq<n; iq++)
		{
			v[ip][iq]=0.0;
		}
		v[ip][ip]=1.0;
	}
	for (ip=0; ip<n; ip++)
	{
		b[ip]=a[ip][ip];
		d[ip]=b[ip];
		z[ip]=0.0;
	}
	nrot=0;
	for (i=0; i<50; i++)
	{
		sm=0.0;
		for (ip=0; ip<n-1; ip++)
		{
			for (iq=ip+1; iq<n; iq++)
			{
				sm=sm+fabs(a[ip][iq]);
			}
		}
		if (sm==0.0)
		{
			return 0;
		}
		if (i<3)
		{
			tresh=0.2*sm/double(n*n);
		}
		else
		{
			tresh=0.0;
		}
		for (ip=0; ip<n-1; ip++)
		{
			for (iq=ip+1; iq<n; iq++)
			{
				g=100.0*fabs(a[ip][iq]);
				sss=fabs(d[ip])+g;
				ddd=fabs(d[iq])+g;
				if ((i>3)&&(sss==fabs(d[ip]))&&(ddd==fabs(d[iq])))
				{
					a[ip][iq]=0.0;
				}
				else
				{
					if (fabs(a[ip][iq])>tresh)
					{
						h=d[iq]-d[ip];
						if ((fabs(h)+g)==fabs(h))
						{
							t=a[ip][iq]/h;
						}
						else
						{
							theta=0.5*h/a[ip][iq];
							t=1.0/(fabs(theta)+sqrt(1.0+pow(theta, 2)));
							if (theta<0.0)
								t=-t;
						}
						c=1.0/sqrt(1.0+t*t);
						s=t*c;
						tau=s/(1.0+c);
						h=t*a[ip][iq];
						z[ip]=z[ip]-h;
						z[iq]=z[iq]+h;
						d[ip]=d[ip]-h;
						d[iq]=d[iq]+h;
						a[ip][iq]=0.0;
						for (j=0; j<=ip-1; j++)
						{
							g=a[j][ip];
							h=a[j][iq];
							a[j][ip]=g-s*(h+g*tau);
							a[j][iq]=h+s*(g-h*tau);
						}
						for (j=ip+1; j<=iq-1; j++)
						{
							g=a[ip][j];
							h=a[j][iq];
							a[ip][j]=g-s*(h+g*tau);
							a[j][iq]=h+s*(g-h*tau);
						}
						for (j=iq+1; j<n; j++)
						{
							g=a[ip][j];
							h=a[iq][j];
							a[ip][j]=g-s*(h+g*tau);
							a[iq][j]=h+s*(g-h*tau);
						}
						for (j=0; j<n; j++)
						{
							g=v[j][ip];
							h=v[j][iq];
							v[j][ip]=g-s*(h+g*tau);
							v[j][iq]=h+s*(g-h*tau);
						}
						nrot=nrot+1;
					}
				}
			}
		}
		for (ip=0; ip<n; ip++)
		{
			b[ip]=b[ip]+z[ip];
			d[ip]=b[ip];
			z[ip]=0.0;
		}
	}
	return 1;
}

int Eigs(double **a, int n, double **v, double eps, int jt)
{
	int i,j,p,q,l;
    double fm,cn,sn,omega,x,y,d;
    l=1;
    for (i=0; i<=n-1; i++)
	{ 
		v[i][i]=1.0;
		for (j=0; j<=n-1; j++)
		{
			if (i!=j) 
			{
				v[i][j]=0.0;
			}
		}
	}
	while (1==1)
	{ 
		fm=0.0;
		for (i=0; i<=n-1; i++)
		{
			for (j=0; j<=n-1; j++)
			{ 
				d=fabs(a[i][j]);
				if ((i!=j)&&(d>fm))
				{ 
					fm=d;
					p=i;
					q=j;
				}
			}
		}
		if (fm<eps)  
		{
			return(1);
		}
		if (l>jt)  
		{
			return(-1);
		}
		l=l+1;
		x=-a[p][q];
		y=(a[q][q]-a[p][p])/2.0;
		omega=x/sqrt(x*x+y*y);
		if (y<0.0)
		{
			omega=-omega;
		}
		sn=1.0+sqrt(1.0-omega*omega);
		sn=omega/sqrt(2.0*sn);
		cn=sqrt(1.0-sn*sn);
		fm=a[p][p];
		a[p][p]=fm*cn*cn+a[q][q]*sn*sn+a[p][q]*omega;
		a[q][q]=fm*sn*sn+a[q][q]*cn*cn-a[p][q]*omega;
		a[p][q]=0.0;
		a[q][p]=0.0;
		for (j=0; j<=n-1; j++)
		{
			if ((j!=p)&&(j!=q))
			{ 
				fm=a[p][j];
				a[p][j]=fm*cn+a[q][j]*sn;
				a[q][j]=-fm*sn+a[q][j]*cn;
			}
		}
		for (i=0; i<=n-1; i++)
		{
			if ((i!=p)&&(i!=q))
			{ 
				fm=a[i][p];
				a[i][p]=fm*cn+a[i][q]*sn;
				a[i][q]=-fm*sn+a[i][q]*cn;
			}
		}
		for (i=0; i<=n-1; i++)
		{ 
			fm=v[i][p];
			v[i][p]=fm*cn+v[i][q]*sn;
			v[i][q]=-fm*sn+v[i][q]*cn;
		}
	}
	return(1);
}