// LU_hinagata.c 
// LU 分解法の雛形プログラム

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

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

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

void matprint(double a[][N]);	// 行列を出力する関数
void vecprint(double *b);		// ベクトルを出力する関数

main()
{
	double a[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;		// 誤差測定のための変数

	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++) b[i] = i + 1;
	// 配列番号が 0 から始まるのでプリントの式とは違っている

	///// 初期設定(ランダムな行列・ベクトル) /////

//	for(i = 0; i < size; i++) for(j = 0; j < size; j++) a[i][j] = REALRANDOM(10);
//	for(i = 0; i < size; i++) b[i] = REALRANDOM(10);

	///// 初期設定(ランダムな整数行列・ベクトル) /////

//	for(i = 0; i < size; i++) for(j = 0; j < size; j++) a[i][j] = RANDOM(10);
//	for(i = 0; i < size; i++) b[i] = RANDOM(10);

	///// 初期設定の出力 /////

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

	///// LU分解 /////

	// 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 分解

//////////////////////////////////////////
////// 未完成部分: L と U を計算せよ /////
//////////////////////////////////////////

	printf("----- LU分解 -----\n\n");
	// L, U を出力
	printf("L:\n"); matprint(ell); printf("\n");
	printf("U:\n"); matprint(u); printf("\n");

	// 検算のため積 LU を計算
	for(i = 0; i < size; i++){
		for(j = 0; j <= size; j++){
			w = 0.0;
			for(k = 0; k < size; k++) w += ell[i][k] * u[k][j];
			c[i][j] = w;
		}
	}
	printf("----- 検算-----\n\n");
	printf("LU:\n"); matprint(c); printf("\n");
	printf("A:\n"); matprint(a); printf("\n");

	///// 前進代入:  Ly=b の解 y の計算 /////

	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];
	}
	printf("----- 前進代入 -----\n\n");
	printf("y:\n"); vecprint(y); printf("\n");

	///// 後退代入:  Ux=y の解 x の計算 /////

	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;
	}
	printf("----- 後退代入 -----\n\n");
	printf("x:\n"); vecprint(x); printf("\n");

	// 検算のため積 Ax を計算
	for(i = 0; i < size; i++){
		w = 0.0;
		for(k = 0; k < size; k++) w += a[i][k] * x[k];
		z[i] = w;
	}
	printf("----- 検算 -----\n\n");
	printf("Ax:\n"); vecprint(z); printf("\n");
	printf("b:\n"); vecprint(b); printf("\n");

	e = 0.0;
	for(i = 0; i < size; i++){
		f = fabs(z[i] - b[i]);
		if(f > e) e = f;
	}
	printf("最大誤差 = %12.10lf\n\n", e);
}


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

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

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