// pm_hinagata.c
// Power Method (累乗法)による最大固有値の計算
// 雛形プログラム
// 高知大学理工学部情報科学ス専門科目　数値解析（塩田）

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

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

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

#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 initialize(double a[][N]);

// 累乗法を用いて行列 A の絶対値最大の固有値と、
// それに対する固有ベクトル、
// 反復回数を返す関数
void maxeigen(double a[][N], double *e, double v[], int *c);

// maxeigen() の出力と、検算を行う関数
void check(double a[][N], double e, double v[], int c);

// ランダムな単位ベクトルを返す関数
void randomunitvec(double u[]);

// 列ベクトルをコピーする関数
void veccopy(double u[], double v[]);

// 行列と列ベクトルの積を計算する関数
void matvecmul(double a[][N], double u[], double v[]);

// 列ベクトルの和を計算する関数
void vecadd(double u[], double v[], double w[]);

// 列ベクトルの差を計算する関数
void vecsub(double u[], double v[], double w[]);

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

// 列ベクトルの長さを返す関数
double norm(double v[]);

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

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

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

// 転置行列を返す関数
void transpose(double a[][N], double b[][N]);

// 列ベクトルのスカラー倍を計算する関数
void vecscmul(double c, double u[], double v[]);

// 行列のスカラー倍を計算する関数
void matscmul(double c, double a[][N], double b[][N]);

// 列ベクトルと行ベクトルの積を計算する関数
void vecmul(double u[], double v[], double a[][N]);

// 行列の差を計算する関数
void matsub(double a[][N], double b[][N], double c[][N]);

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

int main()
{
	double a[N][N];	// 行列 A を表す２次元配列
	double v1[N];	// 終了時に固有ベクトルとなるベクトル v
	double e1;		// 固有値を格納する変数
	int c;			// 反復回数を数えるカウンタ
	clock_t t0, t1;	// 実行時間を測る変数

	initialize(a);

	t0 = clock();	// 開始時刻
	maxeigen(a, &e1, v1, &c);
	t1 = clock();	// 終了時刻
	check(a, e1, v1, c);
	printf("処理時間 = %d ms\n\n", t1 - t0);
	return 0;
}

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


// 行列 A の初期設定
void initialize(double a[][N])
{
	int i, j, mode;
	srand(time(NULL));	// 乱数の初期化

	printf("1: 2次行列の例\n");
	printf("2: ヒルベルト行列\n");
	printf("3: ランダムな整数係数対称行列\n");
	printf("4: ランダムな実数係数対称行列\n");
	printf("5: ランダムな整数係数行列（必ずしも収束しない）\n");
	printf("6: ランダムな実数係数行列（必ずしも収束しない）\n");
	scanf("%d", &mode);
	getchar();

	switch(mode){
		case 1: // 2次行列の例
			size = 2;
			a[0][0] = 3; a[0][1] = 2; a[1][0] = 1; a[1][1] = 2;
			break;

		case 2: // ヒルベルト行列
			printf("行列のサイズは : ");
			scanf("%d", &size);
			getchar();
			for(i = 0; i < size; i++) for(j = 0; j < size; j++) a[i][j] = 1.0 / (i + j + 1);
			break;

		case 3:	// ランダムな整数係数対称行列
			printf("行列のサイズは : ");
			scanf("%d", &size);
			getchar();
			for(i = 0; i < size; i++) for(j = i; j < size; j++){ a[i][j] = RANDOM(10); a[j][i] = a[i][j]; }
			break;

		case 4: // ランダムな実数係数対称行列
			printf("行列のサイズは : ");
			scanf("%d", &size);
			getchar();
			for(i = 0; i < size; i++) for(j = 0; j < size; j++){ a[i][j] = REALRANDOM(10); a[j][i] = a[i][j];}
			break;

		case 5:	// ランダムな整数係数行列
			printf("行列のサイズは : ");
			scanf("%d", &size);
			getchar();
			for(i = 0; i < size; i++) for(j = 0; j < size; j++) a[i][j] = RANDOM(10);
			break;

		case 6: // ランダムな実数係数行列
			printf("行列のサイズは : ");
			scanf("%d", &size);
			getchar();
			for(i = 0; i < size; i++) for(j = 0; j < size; j++) a[i][j] = REALRANDOM(10);
			break;
	}

	printf("A:\n"); matprint(a); printf("\n");
	return;
}

// 累乗法を用いて行列 A の絶対値最大の固有値 e と、
// それに対する固有ベクトル v、
// 反復回数 c を返す関数
void maxeigen(double a[][N], double *e, double v[], int *c)
{
	double u[N];		// 直前のベクトルを格納するベクトル u
	double w[N];		// u + v, u - v を計算するためのベクトル変数
	double pnorm;		// u + v の長さ
	double mnorm;		// u - v の長さ
	int i, j;			// for文の制御変数
	double eps = 1e-8;	// 許容誤差

	///// 固有ベクトルの計算 /////

	// v = ランダムな単位ベクトル
	randomunitvec(v);

	*c = 0;
	do{
		////////////////////////////////////////////
		///// 課題：ここで以下の処理を行え     /////
		/////     　　　　　　　　　　　　　   /////
		///// (1) v を u にコピーする　　　    /////
		///// (2) Au を v に格納する　　　　   /////
		///// (3) v を単位ベクトルに取り替える /////
		////////////////////////////////////////////

		(*c)++;
		if(*c > 256){
			printf("256 回反復しましたが収束しませんでした。\n");
			exit(0);
		}
		if(trace){printf("v:\n"); vecprint(v); printf("\n");}
		// u と v を比較するために
		// u + v, u - v の長さをそれぞれ計算
		vecadd(u, v, w); pnorm = norm(w);
		vecsub(u, v, w); mnorm = norm(w);
	}
	// pnorm または mnorm が許容誤差以下ならば
	// 方向が安定したことになるので反復を終了
	while((pnorm > eps) && (mnorm > eps));

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

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

	return;
}

// maxeigen() の出力と、検算を行う関数
void check(double a[][N], double e, double v[], int c)
{
	double w[N], x[N], y[N];

	///// 計算結果出力 /////

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

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

	printf("----- 検算-----\n\n");
	matvecmul(a, v, w);
	printf("Av:\n"); vecprint(w); printf("\n");
	vecscmul(e, v, x);
	printf("λv:\n"); vecprint(x); printf("\n");
	vecsub(w, x, y);
	printf("Av - λv の長さ = %12.10f\n", norm(y));

	return;
}

// ランダムな単位ベクトルを返す関数
void randomunitvec(double u[])
{
	int i;
	srand(time(NULL));	// 乱数の初期化
	for(i = 0; i < size; i++) u[i] = REALRANDOM(1);
	unitvector(u);
	return;
}

// 列ベクトルをコピーする関数
// v = u;
void veccopy(double u[], double v[])
{
	int i;
	for(i = 0; i < size; i++) v[i] = u[i];
	return;
}

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

// 列ベクトルの和を計算する関数
// w = u + v
void vecadd(double u[], double v[], double w[])
{
	int i;
	for(i = 0; i < size; i++) w[i] = u[i] + v[i];
	return;
}

// 列ベクトルの差を計算する関数
// w = u - v
void vecsub(double u[], double v[], double w[])
{
	int i;
	for(i = 0; i < size; i++) w[i] = u[i] - v[i];
	return;
}

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

// 列ベクトルの長さを返す関数
double norm(double v[])
{
	return sqrt(innerproduct(v, v));
}

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

// 行列を出力する関数
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;
}

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

// 転置行列を返す関数
// B = tA
void transpose(double a[][N], double b[][N])
{
	int i, j;
	for(i = 0; i < size; i++)
		for(j = 0; j < size; j++)
			b[i][j] = a[j][i];
	return;
}

// 列ベクトルのスカラー倍を計算する関数
// v = cu
void vecscmul(double c, double u[], double v[])
{
	int i;
	for(i = 0; i < size; i++) v[i] = c * u[i];
	return;
}

// 行列のスカラー倍を計算する関数
// B = cA
void matscmul(double c, double a[][N], double b[][N])
{
	int i, j;
	for(i = 0; i < size; i++)
		for(j = 0; j < size; j++)
			b[i][j] = c * a[i][j];
	return;
}


// 列ベクトルと行ベクトルの積を計算する関数
// A = uv
void vecmul(double u[], double v[], double a[][N])
{
	int i, j;
	for(i = 0; i < size; i++)
		for(j = 0; j < size; j++)
			a[i][j] = u[i] * v[j];
	return;
}

// 行列の差を計算する関数
// C = A - B
void matsub(double a[][N], double b[][N], double c[][N])

{	int i, j;
	for(i = 0; i < size; i++)
		for(j = 0; j < size; j++)
			c[i][j] = a[i][j] - b[i][j];
	return;
}