// LU_hinagata.c 
// LU 分解法雛形プログラム
//
// 163行目付近の課題その１、175行目付近の課題その２を完成せよ

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

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

#define N 128										// 行列の最大サイズ
#define RANDOM(n) (rand()%(2*n+1)-n)				// 絶対値 |n| 以下の整数乱数 
#define REALRANDOM(x) (x*((2.0*rand()/RAND_MAX)-1))	// 絶対値 |x| 以下の実数乱数 

///// プロトタイプ宣言 /////

// 行列 A, ベクトル b の初期設定
void initialize(double a[][N], double b[]);	
// LU分解法の Step 1（LU分解）
void LUstep1(double a[][N], double ell[][N], double u[][N]);
// LU分解法の Step 2（前進代入）
void LUstep2(double ell[][N], double y[], double b[]);
// LU分解法の Step 3（後退代入）
void LUstep3(double u[][N], double x[], double y[]);
// 行列の積を計算する関数
void matmul(double a[][N], double b[][N], double c[][N]);
// 行列と列ベクトルの積を計算する関数
void matvecmul(double a[][N], double b[], double c[]);
// 行列を出力する関数
void matprint(double a[][N]);
// 列ベクトルを出力する関数
void vecprint(double b[]);
// 行列をコピーする関数
void matcopy(double a[][N], double b[][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;		// 誤差測定のための変数

	// 初期設定の出力
	initialize(a, b);
	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, ベクトル b の初期設定
void initialize(double a[][N], double b[N])
{
	int i, j, mode;
	srand(time(NULL));	// 乱数の初期化

	printf("行列のサイズは？\n");
	scanf("%d", &size);

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

	switch(mode){
		case 1:
			for(i = 0; i < size; i++){
				for(j = 0; j < size; j++)
					a[i][j] = 1.0 / (i + j + 1);
				b[i] = 1.0;
			}
			break;

		case 2:
			for(i = 0; i < size; i++){
				for(j = 0; j < size; j++)
					a[i][j] = REALRANDOM(10);
				b[i] = 1.0;
			}
			break;

		case 3:
			for(i = 0; i < size; i++){
				for(j = 0; j < size; j++)
					a[i][j] = RANDOM(10);
				b[i] = 1;
			}
			break;
	}
	return;
}

// 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;
}

// LU分解法の Step 2（前進代入）
// Ly = b の解 y の計算
void LUstep2(double ell[][N], double y[], double b[])
{
    //////////////////////////////////////////////
    /////     課題その２：ここを完成せよ     /////
    //////////////////////////////////////////////

	return;
}

// LU分解法の 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;
}

// 行列の積を計算する関数
// C = AB
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;
}

// 行列と列ベクトルの積を計算する関数
// c = Ab
void matvecmul(double a[][N], double b[], double c[])
{
	int i, j;
	double w;
	for(i = 0; i < size; i++){
		w = 0.0;
		for(j = 0; j < size; j++) w += a[i][j] * b[j];
		c[i] = w;
	}
	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 b[])
{
	int i;
	for(i = 0; i < size; i++)
		printf("%15.10lf\n", b[i]);
	return;
}

// 行列をコピーする関数
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;
}