// rep10_hinagata.c
// 累乗法による最大固有値の計算
// 雛形プログラム

#include<stdio.h>
#include<stdlib.h>
#include<math.h>

// 行列のサイズは大域変数としてここで宣言する。
int	size = 4;

////////// マクロ等の定義 //////////

#define trace 0	// 途中経過表示 -> 1 / 非表示 -> 0 に設定せよ
#define N 128	// 行列の最大サイズ
#define RANDOM(n) (rand()%(2*n+1)-n)				// 絶対値 |n| 以下の整数乱数
#define REALRANDOM(x) (x*((2.0*rand()/RAND_MAX)-1))	// 絶対値 |x| 以下の実数乱数 

////////// 関数のプロトタイプ宣言 //////////

// 行列 a を出力する関数
void matprint(double a[][N]);

// ベクトル v を出力する関数
void vecprint(double v[]);		

// 行列 a とベクトル v の積をベクトル w に格納する関数
void matvecmul(double a[][N], double v[], double w[]);

// ベクトル u と v の内積を返す関数
double innerproduct(double u[], double v[]);

// ベクトル v を単位ベクトルに取り替えて上書きする関数
void unitvector(double v[]);

////////// メインルーティン //////////

main()
{
	double a[N][N];		// 行列 A を表す２次元配列
	double v[N];		// 終了時に固有ベクトルとなるベクトル v
	double u[N];		// 直前のベクトルを格納するベクトル u
	double pmax;		// u + v の成分の絶対値の最大値
	double mmax;		// u - v の成分の絶対値の最大値
	double w;		// 使い回しする変数
	double eigenvalue;	// 固有値を格納する変数
	int i, j;		// for文の制御変数
	int c = 0;		// 反復回数を数えるカウンタ
	double adm = 1e-8;	// 許容誤差

	srand(time(NULL));	// 乱数の初期化

////////// 初期設定 ///////////

	// 課題の行列の設定
	for(i = 0; i < size; i++) for(j = 0; j < size; j++) a[i][j] = 1.0 / (i + j + 1);
	// ランダムな整数行列の設定
//	for(i = 0; i < size; i++) for(j = 0; j < size; j++) a[i][j] = RANDOM(10);
	// ランダムな行列の設定
//	for(i = 0; i < size; i++) for(j = 0; j < size; j++) a[i][j] = REALRANDOM(10);
	printf("A:\n"); matprint(a); printf("\n");

	// ランダムなベクトル v の初期化
	for(i = 0; i < size; i++) v[i] = REALRANDOM(1);
	// v を単位ベクトルに取り替える
	unitvector(v);

////////// 反復部分 ///////////

	do{
		// v を u にコピー
		for(i = 0; i < size; i++) u[i] = v[i];
		// Au を v に格納
		matvecmul(a, u, v);
		// v を単位ベクトルに取り替える
		unitvector(v);
		// u と v を比較するために
		// u + v, u - v の成分の絶対値の最大値をそれぞれ計算
		pmax = 0.;
		mmax = 0.;
		for(i = 0; i < size; i++){
			w = fabs(v[i] + u[i]);
			if(w > pmax) pmax = w;
			w = fabs(v[i] - u[i]);
			if(w > mmax) mmax = w;
		}
		c++;
		if(trace){printf("v:\n"); vecprint(v); printf("\n");}
		// pmax または mmax が許容誤差以下ならば
		// 方向が安定したことになるので反復を終了
	}
	while((pmax > adm) && (mmax > adm));

////////// 固有値の計算 ///////////

	// もう一度 Av を計算して u に格納
	matvecmul(a, v, u);
	// 固有値は内積 (Av, v) と (v, v) の比として求める

	//////////////////////////////////////////////////////////
	////// 未完成部分: 固有値 eigenvalue の計算式を書け //////
	//////////////////////////////////////////////////////////

////////// 計算結果 ///////////

	printf("----- 計算結果-----\n\n");
	printf("固有値 λ:\n");
	printf("%12.10lf\n", eigenvalue);
	printf("\n");
	printf("固有ベクトル v:\n"); vecprint(v); printf("\n");
	printf("反復回数 = %d\n\n", c);

////////// 検算 ///////////

	printf("----- 検算-----\n\n");
	printf("Av:\n"); vecprint(u); printf("\n");

	// 検算のため v の固有値倍を u に格納
	for(i = 0; i < size; i++)
		u[i] = eigenvalue * v[i];
	printf("λA:\n"); vecprint(u); printf("\n");
}


////////// 関数定義 //////////

// 行列 a を出力する関数
void matprint(double a[][N])
{
	int i, j;
	for(i = 0; i < size; i++){
		for(j = 0; j < size; j++) printf("%12.7lf ", a[i][j]);
		printf("\n");
	}
	return;
}

// ベクトル b を出力する関数
void vecprint(double b[])
{
	int i;
	for(i = 0; i < size; i++) printf("%12.7lf\n", b[i]);
	return;
}

// 行列 a とベクトル v の積をベクトル w に格納する関数
void matvecmul(double a[][N], double v[], double w[])
{
	int i, j;
	double z;
	for(i = 0; i < size; i++){
		z = 0.;
		for(j = 0; j < size; j++) z += a[i][j] * v[j];
		w[i] = z;
	}
	return;
}

// ベクトル u と v の内積を返す関数
double innerproduct(double u[], double v[])
{
	int i;
	double w = 0.;
	for(i = 0; i < size; i++) w += u[i] * v[i];
	return w;
}

// ベクトル v を単位ベクトルに取り替えて上書きする関数
void unitvector(double v[])
{
	int i;
	double len;
	len = sqrt(innerproduct(v, v));	
	for(i = 0; i < size; i++) v[i] /= len;
	return;
}