/* power2.c 
数値解析(C)
累乗法による固有値の計算
gcc power2.c -lm
*/

#include"sqmatrix.h"
// sqmatrix.h を同じディレクトリに download せよ

// 最大固有値 lambda と、lambda に対する固有ベクトル v の計算
void findmaxeigen(MATRIX a, double *lambda, VECTOR *v)
{
	VECTOR v0,v1,w;
	double n,e1,e2,err=1.0e-10;

	// ランダムな単位ベクトル v1 （初期値）の生成 
	v1=randomvec(1.0);
	n=vecnorm(v1);
	v1=vecscmul(1.0/n,v1);

	/* 新 v <- ( Av 方向の単位ベクトル ) の反復 */
	do{
		v0=v1;     // 直前の v1 を v0 として保存 
		v1=matvecmul(a,v0);         // 新v1 = A 旧v
		n=vecnorm(v1);              // v1 の長さを測って
		v1=vecscmul(1.0/n,v1);      // v1 を長さ 1 に正規化
		e1=vecnorm(vecsub(v1,v0));  // v1 と  v0 の差を測定
		e2=vecnorm(vecadd(v1,v0));  // v1 と -v0 の差を測定
	}while(e1>err && e2>err);           // 方向が安定したら終了
	*v = v1;

	/* 固有値の計算 */
	w=matvecmul(a,v1);
	*lambda=innerprod(w,v1)/innerprod(v1,v1);
	return;
}

main()
{
	MATRIX a,a0;
	VECTOR u,v;
	int k;
	double lambda,x;

	RANDOMIZE();  // 乱数の初期化
	SIZE=6;	      // 行列のサイズの設定

/*
	// ランダムな対称行列 A の生成
	a=randommat(5.0);
	a=matadd(a,transmat(a));
*/

/*
	// ランダムな整数係数対称行列 A の生成
	a=randomintmat(5);
	a=matadd(a,transmat(a));
*/

	// 課題の行列の生成

////////// ここを埋めよ //////////

	// 行列 A の表示
	printf("A :\n"); matprint(a); printf("\n");
	// 検算のために A を A0 として保存
	a0=a;

	for(k=0;k<SIZE;k++){
		// 最大固有値と、それに対する固有ベクトルの計算
		findmaxeigen(a,&lambda,&v);
	
		/* 出力と検算 */
		printf("----------------\n\n"); 
		printf("固有値 λ = %lf\n\n",lambda);
		printf("λに対する固有ベクトル v  :\n"); vecprint(v); printf("\n");
		printf("----- 検算 -----\n\n"); 
		printf("Av :\n"); vecprint(matvecmul(a0,v)); printf("\n");
		printf("λv :\n"); vecprint(vecscmul(lambda,v)); printf("\n");

		// デフレーション
		// A の転置行列の最大固有値 x と固有ベクトル u を求める
		findmaxeigen(transmat(a),&x,&u);
		// u と v の内積が 1 になるように u の長さを調整
		u=vecscmul(1.0/innerprod(u,v),u);
		// 新A = 旧A - lambda v (転置 u)
		a=matsub(a,matscmul(lambda,vecrowvecmul(v,transvec(u))));
	}
}

/* 実行例

A :
  -8.0000000   0.0000000   5.0000000
   0.0000000   2.0000000   2.0000000
   5.0000000   2.0000000   4.0000000

----------------

固有値 λ = -9.850128

λに対する固有ベクトル v  :
  -0.9362495
  -0.0584696
   0.3464363

----- 検算 -----

Av :
   9.2221778
   0.5759333
  -3.4124418

λv :
   9.2221778
   0.5759333
  -3.4124418

----------------

固有値 λ = 6.586147

λに対する固有ベクトル v  :
  -0.2997628
  -0.3813558
  -0.8744770

----- 検算 -----

Av :
  -1.9742821
  -2.5116656
  -5.7594337

λv :
  -1.9742821
  -2.5116656
  -5.7594337

----------------

固有値 λ = 1.263981

λに対する固有ベクトル v  :
   0.1832458
  -0.9225774
   0.3395172

----- 検算 -----

Av :
   0.2316193
  -1.1661203
   0.4291433

λv :
   0.2316193
  -1.1661203
   0.4291433

*/