/*
 * ガウスの消去法による連立一次方程式の求解
 * 基本プログラム
 * コンパイル：cc gauss.c -o gauss -lm
 *
 * 問題：連立一次方程式 A x = b
 *   A：係数行列（既知）
 *   b：定数ベクトル（既知）
 *   x：変数ベクトル（未知）
 */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

// 問題番号
#define Q 1

// 計算方法オプション
// #define USE_DOUBLE	// 倍精度実数を利用する
// #define USE_PIVOT	// ピボット選択を利用する
// #define USE_SCALE	// スケーリングを利用する

// 表示方法オプション
#define TRACE		// 途中経過を表示する

// #define COUNT	// （上級者向け）計算量を表示する
// int count = 0;	// （上級者向け）計算量計測用カウンタ
		// 乗除算ごとに count++ し，最後に count を表示せよ．

#ifdef USE_DOUBLE
	#define PRECISION "[倍精度実数]"
	#define FLOAT double		// 8byte，±1.7×10±308，有効数字15桁 
	#define Abs(x) fabs(x)		// 倍精度絶対値関数
#else
	#define PRECISION "[単精度実数]"
	#define FLOAT float		// 4byte，±3.4×10±38，有効数字7桁 
	#define Abs(x) fabsf(x)		// 単精度絶対値関数
#endif
// 注意：絶対値の計算には，fabs() や fabsf() ではなく，この Abs() を使用せよ!!

#ifdef USE_PIVOT
	#define PIVOT "[ピボット選択あり]"
	#ifdef USE_SCALE
		#define SCALE "[スケーリングあり]"
	#else
		#define SCALE "[スケーリングなし]"
	#endif
#else
	#define PIVOT "[ピボット選択なし]"
	#define SCALE ""
#endif

#if Q == 1
	#define PROBLEM "動作テスト用問題(1) {厳密解：x1=-1.0, x2=1.0}"
	// どの計算方法でも正確に計算できていることを確認せよ．
	#define ROW 2	// 拡大行列の行数（方程式の元数）
	#define COL 3	// 拡大行列の列数（元数＋１）
	FLOAT a[ROW][COL] = {	// 拡大行列 [A|b]
		{1.0, 2.0,  1.0},
		{2.0, 1.0, -1.0},
	};

#elif Q == 2
	#define PROBLEM "動作テスト用問題(2) {厳密解：x1=0.5, x2=0.5, x3=0.5}"
	// どの計算方法でも正確に計算できていることを確認せよ．
	#define ROW 3	// 拡大行列の行数（方程式の元数）
	#define COL 4	// 拡大行列の列数（元数＋１）
	FLOAT a[ROW][COL] = {	// 拡大行列 [A|b]
		{2.0, 1.0, 1.0, 2.0},
		{1.0, 2.0, 1.0, 2.0},
		{1.0, 1.0, 2.0, 2.0},
	};

#elif Q == 3
	#define PROBLEM "動作テスト用問題(3) {厳密解：x1=0.5, x2=0.5}"
	// 単精度では，ピボット選択なしだと不正確．
	// 倍精度では，ピボット選択なしでも計算可能．
	// ピボット選択の処理を実装し，単精度でも計算可能とせよ．
	#define ROW 2	// 拡大行列の行数（方程式の元数）
	#define COL 3	// 拡大行列の列数（元数＋１）
	FLOAT a[ROW][COL] = {	// 拡大行列 [A|b]
		{1.0e-10, 1.0, 0.5},
		{1.0,     1.0, 1.0},
	};

#elif Q == 4
	#define PROBLEM "動作テスト用問題(4) {厳密解：x1=0.0, x2=1.0, x3=0.0}"
	// 単精度では，ピボット選択・スケーリングなしだと計算不可能．
	// 倍精度では，ピボット選択・スケーリングなしでも計算可能．
	// ピボット選択・スケーリングの処理を実装し，単精度でも計算可能とせよ．
	#define ROW 3	// 拡大行列の行数（方程式の元数）
	#define COL 4	// 拡大行列の列数（元数＋１）
	FLOAT a[ROW][COL] = {	// 拡大行列 [A|b]
		{9.000000e-15,	1.000000e-15,	5.000000e-15,	1.0e-15},
		{9.000001e+00,	1.000000e+00,	5.000000e+00,	1.0e+00},
		{9.000000e+15,	1.000000e+15,	5.000001e+15,	1.0e+15},
	};

/*
// （上級者向け）独自の問題を創作せよ．
// #elif Q == 5
//	#define PROBLEM "独自問題 {厳密解：...}"
//	#define ROW 3	// 拡大行列の行数（方程式の元数）
//	#define COL 4	// 拡大行列の列数（元数＋１）
//	FLOAT a[ROW][COL] = {	// 拡大行列 [A|b]
//		...
//	};
*/
#endif

// エラー出力
void fatal(char *msg)
{
	fprintf(stderr, "!!! 解析失敗 !!! %s\n", msg);
	exit(1);
}

// 拡大行列を要素表示
void printMat(char *title, FLOAT a[ROW][COL])
{
	int	i, j;
	printf("%s:\n", title);
	for (i = 0; i < ROW; i++) {
		for (j = 0; j < COL; j++) {
			printf("%e\t", a[i][j]);
		}
		printf("\n");
	}
	printf("\n");
}

// 行のスケーリング
void scaleRow(FLOAT a[ROW][COL], int k)
{
// 課題：定義せよ．
// ヒント：第k行をスケーリング（範囲[-1,1]に正規化）．
//   係数行列A の行内の絶対値の最大値を求め，
//   拡大行列[A|b] の行内の全要素を割る．
}

// 行の交換
void swapRow(FLOAT a[ROW][COL], int i1, int i2)
{
// 課題：定義せよ．
// ヒント：第i1行と第i2行とを交換．
//   ただし，i1 == i2 の場合，交換は不要．
}

// ピボット選択
void pivot(FLOAT a[ROW][COL], int k)
{
// 課題：定義せよ．
// ヒント：第k行以降の第k列要素について，
//   絶対値Abs(a[i][k]が最大となる行mを探し，第m行と第k行とを交換する．
//   ただし，USE_SCALE の場合，最大値を求める前に，スケーリングしておくこと．
}

// 前進消去過程（任意行列 → 上三角行列）
void forward(FLOAT a[ROW][COL])
{
	int	i, j, k;
	FLOAT	s;
	for (k = 0; k < ROW-1; k++) {
	/* 高精度化のために...
		ピボット要素 a[k][k] が微小値だと誤差が大きいし，
		もし，ゼロだと解けなくなってしまう．
		そこで，ピボット行を入れ替えよう．
	*/
#ifdef USE_PIVOT
		pivot(a, k);
#endif

		if (a[k][k] == 0.0) fatal("forward(): div/0.");
		for (i = k+1; i < ROW; i++) {
			s = a[i][k]/a[k][k];
			for (j = k; j < COL; j++) {
				a[i][j] -= s*a[k][j];
			}
		}
#ifdef TRACE
		printf("forward %d", k); printMat("",a);
#endif
	}
	// 乗除算を３重ループで計算しているので，計算量はO(n^3)

}

// 後退代入過程（上三角行列 → 対角行列）
void backward(FLOAT a[ROW][COL])
{
	int	i, j, k;
	FLOAT	s;
	for (k = ROW-1; k > 0; k--) {
		if (a[k][k] == 0.0) fatal("backward(): div/0.");
		for (i = k-1; i >= 0; i--) {
			s = a[i][k]/a[k][k];
			for (j = k; j < COL; j++) {
				a[i][j] -= s*a[k][j];
			}
		/* 効率化のために...
			これも O(n^3) だが...実はループjは余計．
			係数行列を実際に対角化する必要はない．
			対角化するつもりの計算を
			b列にだけ適用すれば解を算出できる．
			a[i][COL-1] -= s*a[k][COL-1];
			これで２重ループとなり，計算量はO(n^2)で済む!!
			だが，精度よく対角化できたかどうかも気になるので，
			今回は，ループjも実行している．
		*/
		}
#ifdef TRACE
		printf("backward %d", k); printMat("",a);
#endif
	}
}

// 解の算出（対角行列 → 単位行列）
void finish(FLOAT a[ROW][COL])
{
	int	k;
	FLOAT	s;
	for (k = 0; k < ROW; k++) {
		if (a[k][k] == 0.0) fatal("finish(): div/0.");
		s = a[k][k];
		a[k][k] /= s;		// 要するに a[k][k] = 1.0;
		a[k][COL-1] /= s;	// 解 = b[k] = a[k][COL-1]
	}
}

// ガウスの消去法
void gauss(FLOAT a[ROW][COL])
{
	forward(a);
	backward(a);
	finish(a);
	/* 効率化のために...
		forward() および backward() 内で，
		ピボット要素による除算（a[.][.]/a[k][k]）の重複が多い．
		forward() 内であらかじめ，a[k][k] = 1 となるように，
		ピボット行を正規化（a[k][j] /= a[k][k]）しておけば，
		除算の重複を削減できる．
		（この場合，finish() は不要になる．
		forward() で対角要素 1 の上三角行列化し，
		backward() で単位行列化する．）
	*/
}

int main(int argc, char **argv)
{
	printf("ガウスの消去法\n\n");
	printf("計算方法：%s %s %s\n\n", PRECISION, PIVOT, SCALE);

	printMat(PROBLEM, a);
	gauss(a);
	printMat("結果", a);

	return (0);
}