/* EigenPower.c
累乗法による最大固有値の計算

gcc EigenPower.c -lm
*/

#include"sqmatrix.h"

/* v 方向の単位ベクトルを返す関数 */
VECTOR unitvec(VECTOR v)
{
	return vecscmul(1.0/vecnorm(v),v);
}

/* ランダムな単位ベクトル v を返す関数 */
VECTOR randomunitvec()
{
	VECTOR v=randomvec(1);
	return unitvec(v);
}

void PowerMethodStep1(MATRIX a, double *eval, VECTOR *evec)
{
	VECTOR v,v0;
	double e1,e2,err=1e-8;

	v=randomunitvec();
	do{
		v0=v; 	// 直前の v を v0 として保存 
		v=unitvec(matvecmul(a,v0));
			// 新 v <- ( Av 方向の単位ベクトル )
		e1=vecnorm(vecsub(v,v0));  // v と  v0 の差を測定
		e2=vecnorm(vecadd(v,v0));  // v と -v0 の差を測定
	}while(e1>err && e2>err);          // 方向が安定したら終了
	*evec=v;
	*eval=innerprod(v,matvecmul(a,v))/innerprod(v,v);
}

MATRIX PowerMethodStep2(MATRIX a, double eval, VECTOR evec)
{
	MATRIX b=transmat(a),aa;
	VECTOR u;
	double y;
	PowerMethodStep1(b,&y,&u);
	u=vecscmul(1.0/innerprod(u,evec),u);
	aa=matsub(a,matscmul(eval,vecrowvecmul(evec,transvec(u))));
	return aa;
}

void PowerMethod(MATRIX a0, double *eval, VECTOR *evec)
{
	MATRIX a=a0;
	VECTOR v;
	double k;
	int i;
	for(i=0;i<SIZE;i++){
		PowerMethodStep1(a,&k,&v);
		eval[i]=k;
		evec[i]=v;
		a=PowerMethodStep2(a,k,v);
	}
}

main()
{
	MATRIX a;
	VECTOR v[MAXSIZE];
	double k[MAXSIZE];
	int i,j;

	srand(time(NULL));

	SIZE=8;
	a=randommat(5);
	a=matadd(a,transmat(a));
		// ランダムな対称行列
	printf("A :\n"); matprint(a); printf("\n");

	PowerMethod(a,k,v);
	
	/* 出力と検算 */
	for(i=0;i<SIZE;i++){
		printf("--------------------\n\n");
		printf("%d-th Eigenvalue λ = %lf\n\n",i,k[i]);
		printf("Eigenvector v :\n");
		vecprint(v[i]); printf("\n");
/*
		printf("Check :\n\n"); 
		printf("Av :\n"); 
		vecprint(matvecmul(a,v[i])); printf("\n");
		printf("λv :\n"); 
		vecprint(vecscmul(k[i],v[i])); printf("\n");
*/
		printf("Hit Return Key\n");
		getchar();
	}
}

/*　実行例



*/