// GaussianElimination.c
//ガウスの消去法

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

#define trace 1
// 途中経過をトレースするときは 1 に、しないときは 0 に設定せよ

#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], int m, int n);	// 行列を出力する関数
void vecprint(double *v, int n);				// ベクトルを出力する関数

main()
{
	double a0[N][N], b0[N];	// 連立一次方程式の係数
	double a[N][N];			// 拡大係数行列
	double x[N];			// 未知数ベクトル
	double z[N];			// 検算用ベクトル
	int size;				// 行列のサイズ
	int i, j, n, k, p;
	double m, max, t, c, e, f, w;
	char r;

	srand(time(NULL));	// 乱数の初期化

	///// 方程式の設定(鶴亀算の場合) /////

/*
	size = 2; 	// 行列のサイズを設定

	// 乱数で係数を設定
	a0[0][0] =  1.0;
	a0[0][1] =  1.0;
	a0[1][0] =  2.0;
	a0[1][1] =  4.0;
	b0[0]    =  7.0;
	b0[1]    = 22.0;
*/

	///// 方程式の設定(ランダムな整数係数の場合) /////

	size = 3; 	// 行列のサイズを設定

	// 乱数で係数を設定
	for(i = 0; i < size; i++) for(j = 0; j < size; j++)	a0[i][j] = RANDOM(10);
	for(i = 0; i < size; i++) b0[i] = RANDOM(10);

	// ピボット選択の効果を見るときは次の行のコメントアウトをはずす
	//a0[0][0] = 0.000001;

	///// 方程式の設定(ランダムな実数係数の場合) /////

/*
	size = 3; 	// 行列のサイズを設定

	// 乱数で係数を設定
	for(i = 0; i < size; i++) for(j = 0; j < size; j++) a0[i][j] = REALRANDOM(10);
	for(i = 0; i < size; i++) b0[i] = REALRANDOM(10);

	// ピボット選択の効果を見るときは次の行のコメントアウトをはずす
	a0[0][0] = 0.000001;
*/

////////// シンプルなガウスの消去法 //////////

	///// 拡大係数行列の設定 /////
	for(i = 0; i < size; i++){
		for(j = 0; j < size; j++) a[i][j] = a0[i][j];
		a[i][size] = b0[i];
	}

	printf("----- ピボット選択無し -----\n\n");
	printf("拡大系数行列 :\n"); matprint(a, size, size + 1); printf("\n");

	///// 前進消去 /////

	for(k = 0; k < size - 1; k++){
		// i 行目から, k 行目の (a[i][k] / a[k][k])倍を引く
		for(i = k + 1; i < size; i++){
			c = a[i][k] / a[k][k];
			for(j = k; j <= size; j++) a[i][j] -= c * a[k][j];
		}
		if(trace){
			printf("%d 行目を用いた掃き出し\n",k);
			printf("->\n");
			matprint(a, size, size + 1); printf("\n");
			printf("Hit Return Key\n");
			r = getchar();
		}
	}

	///// 後退代入 /////

	for(i = size - 1; i >= 0; i--){
		w = 0.0;
		for(j = i + 1; j < size; j++) w += a[i][j] * x[j];
		x[i] = (a[i][size] - w) / a[i][i];
	}

	///// 検算: z = Ax /////

	for(i = 0; i < size; i++){
		w = 0.0;
		for(j = 0; j < size; j++) w += a0[i][j] * x[j];
		z[i] = w;
	}

	///// 結果出力 /////

	printf("->\n");
	printf("A' :\n"); matprint(a, size, size + 1); 		
	printf("x :\n"); vecprint(x, size); printf("\n");
	printf("検算 Ax :\n"); vecprint(z, size); printf("\n");
	printf("b :\n"); vecprint(b0, size); printf("\n");

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

////////// 部分ピボット選択を用いたガウスの消去法 //////////

	///// 拡大係数行列の設定 /////
	for(i = 0; i < size; i++){
		for(j = 0; j < size; j++) a[i][j] = a0[i][j];
		a[i][size] = b0[i];
	}

	printf("----- ピボット選択あり -----\n\n");
	printf("拡大系数行列 :\n"); matprint(a, size, size + 1); printf("\n");

	///// 前進消去 /////

	for(k = 0; k < size - 1; k++){

		// 部分ピボット選択
		// a[k][k],...,a[size-1][k] の中で
		// 絶対値最大の要素が a[p][k] であるとする
		p = k; max = fabs(a[k][k]);
		for(i = k + 1; i < size; i++){
			m = fabs(a[i][k]);
			if(m > max){ p = i; max = m; }
		}

		// p != k のときは k 行目と p 行目を入れ替える
		if(p != k){
			for(j = k; j <= size; j++){	t = a[k][j]; a[k][j] = a[p][j]; a[p][j] = t; }
			if(trace){
				printf("%d 行目と %d 行目を交換\n",k,p);
				matprint(a, size, size + 1); printf("\n");
				printf("Hit Return Key\n");
				r = getchar();
			}
		}

		// 掃き出し
		for(i = k + 1; i < size; i++){
			// i 行目から, k 行目の (a[i][k] / a[k][k])倍を引く
			c = a[i][k] / a[k][k];
			for(j = k; j <= size; j++) a[i][j] -= c * a[k][j];
		}
		if(trace){
			printf("%d 行目を用いた掃き出し\n",k);
			printf("->\n");
			matprint(a, size, size + 1); printf("\n");
			printf("Hit Return Key\n");
			r = getchar();
		}
	}

	///// 後退代入 /////

	for(i = size - 1; i >= 0; i--){
		w = 0.0;
		for(j = i + 1; j < size; j++) w += a[i][j] * x[j];
		x[i] = (a[i][size] - w) / a[i][i];;
	}

	///// 検算: z = Ax /////

	for(i = 0; i < size; i++){
		w = 0.0;
		for(j = 0; j < size; j++) w += a0[i][j] * x[j];
		z[i] = w;
	}

	///// 結果出力 /////

	printf("->\n");
	printf("A' :\n"); matprint(a, size, size + 1); 		
	printf("x :\n"); vecprint(x, size); printf("\n");
	printf("検算 Ax :\n"); vecprint(z, size); printf("\n");
	printf("b :\n"); vecprint(b0, size); printf("\n");

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

void matprint(double (*a)[N], int m, int n)
{
	int i, j;
	for(i = 0; i < m; i++){
		for(j = 0; j < n; j++) printf("%12.7lf",a[i][j]);
		printf("\n");
	}
}

void vecprint(double *v, int n)
{
	int i;
	for(i = 0; i < n; i++) printf("%12.7lf\n", v[i]);
}