//////////////////////////////////////////////////////////
//
// rep06_hinagata.c  数値解析 2011.1.26 の課題
// 累乗法による最大固有値の計算　雛形プログラム
// 未完成部分１・２を完成せよ
//
//////////////////////////////////////////////////////////

#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#define N 128		// 行列の最大サイズ

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

// 絶対値 |x| 以下の実数乱数 
double myrandom(double x)
{
	return 2 * x * ((1.0 * rand() / RAND_MAX) - 0.5);
}

// 絶対値 |n| 以下の整数乱数 
int myrandominteger(int n)
{
	return rand() % (2 * n) - n;
}

// 行列 A と列ベクトル u の積 Au を v として返す関数
void matrix_vector_mul(double a[][N], double *u, double *v)
{
////////////////////////////////////////////
////// 未完成部分１: v = Au を計算せよ /////
////////////////////////////////////////////
}

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

// ベクトル v を単位ベクトルに取り替える関数
void unit_vector(double *v)
{
	int i;
	double len;
	len = sqrt(inner_product(v, v));	
	for(i = 0; i < size; i++) v[i] /= len;
}

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

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

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 AllowableError = 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);

	// ベクトル v の初期値の設定（ランダム）
	for(i = 0; i < size; i++) v[i] = myrandom(1);

	// v を単位ベクトルに取り替える
	unit_vector(v);

	// ランダムな整数行列の設定
//	for(i=0; i<size; i++) for(j=0; j<size; j++) a[i][j] = myrandominteger(10);
	// ランダムな行列の設定
//	for(i=0; i<size; i++) for(j=0; j<size; j++) a[i][j] = myrandom(10);

	printf("A:\n");
	print_matrix(a);
	printf("\n");

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

	do{

		// v を u にコピー
		for(i = 0; i < size; i++) u[i] = v[i];
		// Au を v に格納
		matrix_vector_mul(a, u, v);
		// v を単位ベクトルに取り替える
		unit_vector(v);

		// u と v を比較するために
		// u + v, u - v の成分の絶対値の最大値をそれぞれ計算
		pmax = 0.0;
		mmax = 0.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++;
		// pmax または mmax が許容誤差以下ならば
		// 方向が安定したことになるので反復を終了
	}
	while(pmax > AllowableError && mmax  > AllowableError);

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

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

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

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

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

//////////////////////////////////////
////// 未完成部分２: 検算をせよ //////
//////////////////////////////////////
}