repowmtd.c

来自「用幂法求矩阵的最大特证值」· C语言 代码 · 共 59 行

#include "math.h"
#include "ordGauss.c"
#include "head.c"

double rePowMethod(double m[C_M][C_N],double *u,double *y,double para)
{
	int i,j,k,flag;
	double b;
	double temp,tag;
	
	u[0]=1;
	for (i=1;i<Matrix_N ;i++ )	u[i]=0;
	k=0;
	flag=0;
	b=0;
	for (i=0;i<Matrix_N ;i++ )	m[convert_i(i,i)][i]-=para;
	do
	{	
		k++;
		temp=0;
		for (i=0;i<Matrix_N ;i++)	temp+=u[i]*u[i];
		temp=sqrt(temp);
		for (i=0;i<Matrix_N ;i++)	y[i]=u[i]/temp;
		flag=ordGauss(m,y,u);
		if (flag==0) exit("0");
		temp=b;
		b=0;
		for (i=0;i<Matrix_N ;i++)	b+=y[i]*u[i];
		tag=fabs(b-temp)/fabs(b);
	}
	while (tag>E);
	return 1/b;
}

/* for test	*/
/*main()
{
	double test[C_M][C_N];
	double u[Matrix_N],y[Matrix_N];
	double b;
	double temp1,temp2;
	int i,j;

	createMatrix_1(&test);
	b=1;
	b=rePowMethod(&test,u,y,0);
	for (i=0;i<Matrix_N ;i++ )
	{	
		temp1=0;temp2=0;
		for (j=i-2;j<=i+2;j++ )
		{
			if (j>=0 && j<Matrix_N)
			{
				temp1+=test[convert_i(i,j)][j]*y[j];
			}
		}
		printf("%-12f,%-12f\n",temp1,b*y[i]);
	}
}*/