// ガウスの消去法で連立一次方程式を解くプログラム

// 4x + y + z = 9
// x + 3y + z = 10
// 2x + y + 5z = 19

#include <stdio.h>

#define SIZE 3

double A[SIZE][SIZE] = {
	{4, 1, 1},
	{1, 3, 1},
	{2, 1, 5}
};

double b[SIZE]={9,10,19};

double c[SIZE];
double x[SIZE];


int main(int argc, char* argv[]){
	int i,j,k;
	int n = SIZE;


	// 外積方式ガウス法

	// LU分解
	for (k=0; k<n; k++) {
		// kで回す行方向ループ

		// 行列Lの処理
		double div_tmp = 1.0 / A[k][k];
		for (i=k+1; i<n; i++) {
			// iで回す行方向ループ

			// l_ik = a_ik / a_kk
			A[i][k] = A[i][k] * div_tmp;   
		}

		// 行列Uの処理
		for (j=k+1; j<n; j++) {
			// jで回す行方向のループ

			// Tips: ここで値を変数に入れたほうが
			// ループ内で参照するよりメモリアクセス回数が減ってお得
			double a_jk = A[j][k];

			for (i=k+1; i<n; i++) {
				// iで回す列方向のループ
				A[j][i] = A[j][i] - A[k][i] * a_jk; 
			}
		}
	}
	// LU分解終了


	// LU分解後のLUを出力
	printf("LU:\n");
	for(i=0;i<n;i++){
		for(j=0;j<n;j++){
			printf("% 2.2f ",A[i][j]);
		}
		printf("\n");
	}
	printf("\n");


	// 前進代入
	// Lc = b
	for(k=0; k<n; k++){
		// 行方向のループ
		c[k] = b[k];
		for(i=0; i<k; i++){
			// 列方向のループ
			c[k] -= A[k][i] * c[i];
		}
	}


	// 後進代入
	// Ux = c
	// k = n - 1
	x[n-1] = c[n-1] / A[n-1][n-1];
	// k = (n - 2) ~ (0)
	for(k = n - 2; k >= 0; k--){
		x[k] = c[k];
		for(i=k+1; i<n; i++){
			// 列方向のループ
			x[k] -= A[k][i] * x[i];
		}
		x[k] = x[k] / A[k][k];
	}

	// xを出力
	printf("x:\n");
	for(i=0;i<n;i++){
		printf("% 2.2f ",x[i]);
	}
	printf("\n");

	return 0;
}