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

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

main()
{
	MATRIX a;
	VECTOR v,v0,w;
	int i,j;
	double n,e1,e2,err=1.0e-10;
	double lambda;

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

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

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

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

	// 行列 A の表示
	printf("A :\n"); matprint(a); printf("\n");

	/* ランダムな単位ベクトル v （初期値）の生成 */
	v=randomvec(1.0);
	n=vecnorm(v);
	v=vecscmul(1.0/n,v);
	
	/* 新 v <- ( Av 方向の単位ベクトル ) の反復 */
	do{
		v0=v; // 直前の v を v0 として保存 
		v=matvecmul(a,v0);   // 新v = A 旧v
		n=vecnorm(v);        // v の長さを測って
		v=vecscmul(1.0/n,v); // v を長さ 1 に正規化
		// 途中経過を出力する場合は次のコメントアウトをはずす
		//printf("v :\n"); vecprint(v); printf("\n");
		e1=vecnorm(vecsub(v,v0));  // v と  v0 の差を測定
		e2=vecnorm(vecadd(v,v0));  // v と -v0 の差を測定
	}while(e1>err && e2>err);      // 方向が安定したら終了

	/* 固有値の計算 */
	w=matvecmul(a,v);
	lambda=innerprod(w,v)/innerprod(v,v);
	
	/* 出力と検算 */
	printf("最大固有値 λ = %lf\n\n",lambda);
	printf("λに対する固有ベクトル v  :\n"); vecprint(v); printf("\n");
	printf("----- 検算 -----\n\n"); 
	printf("Av :\n"); vecprint(w); printf("\n");
	printf("λv :\n"); vecprint(vecscmul(lambda,v)); printf("\n");
}

/*　実行例

A :
   2.0000000   2.0000000   2.0000000
   2.0000000   4.0000000  -5.0000000
   2.0000000  -5.0000000  -2.0000000

最大固有値 λ = 6.962522

λに対する固有ベクトル v  :
  -0.1715420
  -0.8761408
   0.4505005

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

Av :
  -1.1943646
  -6.1001494
   3.1366192

λv :
  -1.1943646
  -6.1001494
   3.1366192

*/