powermethod.c

C语言幂法程序,幂法的主程序,要用到ElemMatOper.c文件中的函数

/* This program implements the power method of finding the dominant
 * eigenvector and eigenvalue of a matrix. (see the description below)
 * 
 * This program is compiled using gcc 3.4.2 with mingw under Windows
 * operating system or gcc 4.* under Ubuntu Linux operating system.
 * It is convenient to use the software gnuplot to visualize the computation 
 * results with the following commands in Linux shell(bash):
 * #building program
 * gcc -c powermethod.c -o powermethod
 * #executing the program
 * ./powermethod | tee dat
 * #visualizing the result
 * gnuplot
 * #in gnuplot
 * plot "dat" with linespoints
 * #end
 *
 * This program was written by Xiaoke Yang @ School of Automation Science and
 * Electrical Engineering, Beihang University. You can copy, modify or 
 * redistribute it freely with this header kept intact. Welcome reports of 
 * bugs, suggestions, etc. [e-mail]
 *
 * Last modified by Xiaoke Yang, Oct.31,2007
 */
#include "ElemMatOper.c"

#undef IMAX				/*IMAX is the maximus iteration count defined in*/
#define IMAX 20000		/*the header file, and can be redefined here*/
#undef N				/*N is the matrix dimension used here, redefined*/
#define N 3
#undef DEBUG			/*DEBUG is a flag of printing the debug information*/
#define DEBUG 1			/*write 0 here to turn off the debug information*/



/* Function: double PowerIter(double A[N][N], double u[N], double dv)
 * Description: This procedure uses the Power method or Power Iteration
 * to find the dominant eigenvalue and eigenvector of a sqare real matrix A.
 * This method can be used when 
 * (1). A has N linearly independent eigenvectors.
 * (2). A has an eigenvalue that is strictly greater in magnitude than its 
 * other eigenvalues.
 * The power iteration algorithm starts with a vector u0 which has a nonzero
 * component in the direction of the eigenvector associated with the dominant
 * eigenvalue and iterates by u(k+1)=A*u(k)/norm(A*u(k)), then a subsequence
 * of u(k) converges to an eigenvector associated to the dominant eigenvalue.
 * Anyway, there are some cases in which this method might fail.
 * Input:	square matrix A[N][N]
 *			initial vector u0[N]
 *			eigenvector dv[N] of A associated with the dominant eigenvalue
 * Output:	dominant eigenvalue of A
 * Calling Syntax: PowerIter(A,u0,dv)
 * */
double PowerIter(double A[N][N],double u0[N],double dv[N]){
	double u[N],y[N]={0,0,0};
	double u1[N],y1[N]={0,0,0};
	double un,de,erre,errv,de1=0;		/*un is the norm of vector u*/
	int i;								/*de is the dominant eigenvalue*/
	veccopy(u0,u);
	/*iteration starts*/
	for(i=0;i<IMAX;i++){
		/*save data in the last iteration for checking stop conditions*/
		veccopy(y,y1);			/*save eigenvector in the last interation*/
		de1=de;					/*save eigenvalue in the last iteration*/
		/*end saving, begin computing the new data*/
		un=vecnorm2(u);			/*get the norm of u*/
		veccopy(y1,y);
		vecmultcon(u,1.0/un,y);	/*calculate y*/
		matmultvec(A,y,u);		/*u=A*y*/
		de=innerprod(y,u);		/*calculate the eigenvalue*/
		vecmultcon(y,de,u1);	/*u1 is used to check whether de*y=A*y or not*/
		/*end computing the new data, checking the stop conditions*/
		/*relative error between the eigenvalues in present and previous iterations*/
		erre=fabs((de-de1)/de);	
		/*relative error between eigenvectors in present and previous iterations*/
		errv=vecdist(y,y1)/vecnorm2(y);
		if(erre<eps){
			printf("\n#Stopping condition with the eigenvalue achieved at iteration %d with error=%e<%e!\n",i+1,erre,eps);
			printf("#The dominant eigenvalue is: lambda=\n");
			printf("%e\n",de);
			printf("#The corresponding eigenvector is: y(k-1)=\n");
			vecdisp(y);
			printf("#u(k)=A*y(k-1)=\n");
			vecdisp(u);
			printf("#lambda*y(k-1)=\n");
			vecdisp(u1);
			printf("#Error with eigenvectors is:\n");
			printf("%e\n",errv);
			printf("#Error with eivenvalues is:\n");
			printf("%e\n",erre);
			break;
		}
		if(DEBUG==1){
			printf("\n#The dominant eigenvalue after iteration %d is: lambda=\n",i+1);
			printf("%e\n",de);
			printf("#The corresponding eigenvector is: y(k-1)=\n");
			vecdisp(y);
			printf("#u(k)=A*y(k-1)=\n");
			vecdisp(u);
			printf("#lambda*y(k-1)=\n");
			vecdisp(u1);
			printf("#Error with eigenvectors is:\n");
			printf("%e\n",errv);
			printf("#Error with eivenvalues is:\n");
			printf("%e\n",erre);
		}
	}
	veccopy(u,dv);				/*return the eigenvector*/
	return de;					/*return the eigenvalue*/
}


int main(){
/*	double A[N][N]={{30,2,1},{0,30.1,4},{0,0,30.11}};
	double u0[N]={10,0.5,10};
	double dv[N];
*/	
	/*the above is a case in which lambda1~lambda2~lambda3*/
	double A[N][N]={{1,2,3},{0,-3,4},{0,0,3}};
	double u0[N]={1,0.5,1};
	double dv[N];
	PowerIter(A,u0,dv);
}