//////////////////////////////////////////////////////////
//
// rep07_hinagata.c  数値解析 2013.1.9 の課題
// LU 分解法 雛形プログラム
// 未完成部分１・２を完成せよ
//
//////////////////////////////////////////////////////////

#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#define N 128

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

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

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

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文の制御変数

////////// 初期設定 ///////////

	// 課題の行列・ベクトルの設定
	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;

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


////////// Step 1: 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  を出力
	printf("L:\n");
	print_matrix(ell);
	printf("\n");

	// U  を出力
	printf("U:\n");
	print_matrix(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");
	print_matrix(c);
	printf("\n");
	printf("A:\n");
	print_matrix(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");
	print_vector(y);
	printf("\n");

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

///////////////////////////////////////
////// 未完成部分２: x を計算せよ /////
///////////////////////////////////////

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