// LU_hinagata.c 
// 連立一次方程式のLU 分解法
// 雛形プログラム

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

////////// マクロと大域変数の宣言 //////////////////////////////////////////////

#define N 128	// 行列のサイズの上限
int	size;		// 行列のサイズ

////////// LU分解法の関数定義 //////////////////////////////////////////////////

// Step 1（LU分解）: A = LU の形に分解する
void LUstep1(double a[][N], double ell[][N], double u[][N])
{
	int i, j;

	// L, U の初期値は単位行列
	for(i = 0; i < size; i++){
		for(j = 0; j < size; j++){
			ell[i][j] = 0;
			u[i][j] = 0;
		}
		ell[i][i] = 1;
		u[i][i] = 1;
	}

	// LU 分解

    ///////////////////////////////////////////////////
    ///// 課題：                                  /////
	///// ここを完成せよ（変数は適宜追加宣言せよ) /////
    ///////////////////////////////////////////////////

	return;
}

// Step 2（前進代入）: Ly = b の解 y の計算
void LUstep2(double ell[][N], double y[], double b[])
{
	int i, j;
	double w;

	for(i = 0; i < size; i++){
		w = 0.0;
		for(j = 0; j < i; j++) w += ell[i][j] * y[j];
		y[i] = (b[i]- w) / ell[i][i];
	}
	return;
}

// Step 3（後退代入）: Ux = y の解 x の計算
void LUstep3(double u[][N], double x[], double y[])
{
	int i, j;
	double w;
	for(i = size - 1; i >= 0; i--){
		w = 0.0;
		for(j = i + 1; j < size; j++) w += u[i][j] * x[j];
		x[i] = y[i] - w;
	}
	return;
}

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

// 行列 a を出力する関数
void matprint(double a[][N]);
// ベクトル v を出力する関数
void vecprint(double v[]);
// 行列 a を行列 b にコピーする関数
void matcopy(double a[][N], double b[][N]);
// 行列 a とベクトル v の積をベクトル w に格納する関数
void matvecmul(double a[][N], double v[], double w[]);
// 行列 a と行列 b の積を行列 c に格納する関数
void matmul(double a[][N], double b[][N], double c[][N]);

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

int main()
{
	double a[N][N];     // 行列 A
	double a0[N][N];    // 行列 A のコピー（検算用）
	double ell[N][N];   // 行列 L
	double u[N][N];     // 行列 U
	double c[N][N];     // 検算のための行列
	double b[N];        // ベクトル b
	double x[N];        // 解ベクトル x
	double y[N];        // 前進代入で求めるベクトル y
	double z[N];        // 検算のためのベクトル
	double w;           // 総和を求めるための変数
	int i, j, k;        // for文の制御変数
	double e, f;		// 誤差測定のための変数
	int mode;

	////////// 方程式の設定 //////////

	printf("行列を設定してください：\n");
	printf("1 サンプルデータ\n");
	printf("2 ランダムな行列\n");
	printf("3 課題のヒルベルト行列\n");
	scanf("%d", &mode);

	switch(mode){
		case 1:
			// サンプルデータ
			size = 2;
			a[0][0] = 3; a[0][1] = -2; a[1][0] = 1; a[1][1] = 3;
			b[0] = 1; b[1] = 4;
			break;

		case 2:
			// ランダムな行列
			srand(time(NULL));	// 乱数の初期化
			printf("ランダムな行列を生成します。\n");
			printf("行列のサイズを入力してください。\n");
			scanf("%d", &size);
			for(i = 0; i < size; i++){
				for(j = 0; j < size; j++) a[i][j] = 20 * ((1.0 * rand()) / RAND_MAX - 0.5);
				b[i] = 1.0;
			}
			break;

		case 3:
			// 課題のヒルベルト行列
			printf("課題の方程式で計算します。\n");
			printf("行列のサイズを入力してください。\n");
			scanf("%d", &size);
			for(i = 0; i < size; i++){
				for(j = 0; j < size; j++) a[i][j] = 1.0 / (i + j + 1);
				b[i] = 1.;
			}
			break;
	}

	printf("\n----- A と b の設定 -----\n\n");
	printf("A:\n"); matprint(a); printf("\n");
	printf("b:\n"); vecprint(b); printf("\n");

	////////// 実行部 //////////

	// 検算のため A を A0 にコピー
	matcopy(a, a0);
	// LU分解
	LUstep1(a, ell, u);
	// 検算のため積 C = LU を計算
	matmul(ell, u, c);
	// 前進代入:  Ly = b の解 y の計算 
	LUstep2(ell, y, b);
	// 後退代入:  Ux = y の解 x の計算
	LUstep3(u, x, y);
	// 検算のため積 z = Ax を計算
	matvecmul(a0, x, z);
	// 最大誤差の計算
	e = 0.0;
	for(i = 0; i < size; i++){
		f = fabs(z[i] - b[i]);
		if(f > e) e = f;
	}
	// 出力
	printf("----- LU分解 -----\n\n");
	printf("L:\n"); matprint(ell); printf("\n");
	printf("U:\n"); matprint(u); printf("\n");
	printf("----- 検算-----\n\n");
	printf("LU:\n"); matprint(c); printf("\n");
	printf("A:\n"); matprint(a); printf("\n");
	printf("----- 前進代入 -----\n\n");
	printf("y:\n"); vecprint(y); printf("\n");
	printf("----- 後退代入 -----\n\n");
	printf("x:\n"); vecprint(x); printf("\n");
	printf("----- 検算 -----\n\n");
	printf("Ax:\n"); vecprint(z); printf("\n");
	printf("b:\n"); vecprint(b); printf("\n");
	printf("最大誤差 = %12.10lf\n\n", e);

	return 0;
}

////////// 関数定義（以下は触らなくてよい） ////////////////////////////////////

// 行列 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 を行列 b にコピーする関数
void matcopy(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[i][j];
	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;
}

// 行列 a と行列 b の積を行列 c に格納する関数
void matmul(double a[][N], double b[][N], double c[][N])
{
	int i, j, k;
	double x;
	for(i = 0; i < size; i++){
		for(j = 0; j < size; j++){
			x = 0.0;
			for(k = 0; k < size; k++) x += a[i][k] * b[k][j];
			c[i][j] = x;
		}
	}
	return;
}