// JGS_hinagata.c
// 連立一次方程式の反復解法: 
// ヤコビ法, ガウス-ザイデル法
// 雛形プログラム

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

// 行列のサイズは大域変数としてここで宣言する。
int	size;

////////// マクロ等の定義 //////////

#define trace 0	// 途中経過表示 -> 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| 以下の実数乱数 

////////// 関数のプロトタイプ宣言 //////////

// 行列 a を出力する関数
void matprint(double a[][N]);

// ベクトル v を出力する関数
void vecprint(double v[]);	
	
// 行列 a とベクトル v の積をベクトル w に格納する関数
void matvecmul(double a[][N], double v[], double w[]);

// ベクトル u とベクトル v の差 u - v をベクトル w に格納する関数
void vecsub(double u[], double v[], double w[]);

// ベクトル v の長さを返す関数
double norm(double v[]);

////////// ヤコビ法･ガウス-ザイデル法の関数宣言 //////////

// ヤコビ法の１ステップ
// 近似ベクトル x を更新する関数
void Jacobi(double a[][N], double b[], double x[])
{

//////////////////////////////////////////////////////
////// 未完成部分：                              /////
////// 次のガウス-ザイデル法をまねてここを埋めよ /////
//////////////////////////////////////////////////////

}

// ガウス-ザイデル法の１ステップ
// 近似ベクトル x を更新する関数
void GaussSeidel(double a[][N], double b[], double x[])
{
	int i, j;
	double w[N]; 	// w = 新x
	double z;
	for(i = 0; i < size; i++){
		z = b[i];
		for(j = 0; j < i; j++) z -= a[i][j] * w[j];
		for(j = i + 1; j < size; j++) z -= a[i][j] * x[j];
		w[i] = z / a[i][i];
	}
	for(i = 0; i < size; i++) x[i] = w[i];	// w を x に上書き
	return;
}

////////// メインルーティン //////////

main()
{
	double a[N][N];		// 行列 A を表す２次元配列
	double b[N];		// ベクトル b を表す配列
	double x[N];		// 解ベクトル x を表す配列
	double w[N];		// ベクトル Ax を格納する配列
	double c[N];		// 誤差ベクトル Ax - b を格納する配列
	double e;			// 誤差ベクトル Ax - b の長さを格納する変数
	int maxrec = 128;	// 最大反復回数
	double adm = 1e-7;	// 許容誤差
	int i, j, step;

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

	////////// 方程式の設定 //////////

	switch(5){	// 番号を選べ
		default: // 収束する例
			size = 2;
			a[0][0] = 3; a[0][1] = -2; a[1][0] = 1; a[1][1] = 3;
			b[0] = 1; b[1] = 4;
			break;
		case 1:	// 収束しない例
			size = 2;
			a[0][0] = 1; a[0][1] = 3; a[1][0] = 3; a[1][1] = -2;
			b[0] = 4; b[1] = 1;
			break;
		case 2:	// 式・未知数を入れ替えることにより収束の速さが変わる例の１
			size = 3;
			a[0][0] = 3; a[0][1] = 2; a[0][2] =  2;
			a[1][0] = 2; a[1][1] = 5; a[1][2] = -4;
			a[2][0] = 1; a[2][1] = 1; a[2][2] =  7;
			b[0] = 7; b[1] = 3; b[2] = 9;
			break;
		case 3 : // 式・未知数を入れ替えることにより収束の速さが変わる例の２
			size = 3;
			a[0][0] =  7; a[0][1] = 1; a[0][2] = 1;
			a[1][0] = -4; a[1][1] = 5; a[1][2] = 2;
			a[2][0] =  2; a[2][1] = 2; a[2][2] = 3;
			b[0] = 9; b[1] = 3; b[2] = 7;
			break;
		case 4: // ランダムな整数係数の場合
			size = 4;
			for(i = 0; i < size; i++){
				for(j = 0; j < size; j++) a[i][j] = RANDOM(10);
				a[i][i] = RANDOM(10 * size);	// 対角成分は大きく
				b[i] = RANDOM(10);
			}
			break;
		case 5:	// ランダムな実数係数の場合
			size = 4;
			for(i = 0; i < size; i++){
				for(j = 0; j < size; j++) a[i][j] = REALRANDOM(10);
				a[i][i] *= size;	// 対角成分は大きく
				b[i] = REALRANDOM(10);
			}
			break;
		case 6:	// 課題の方程式
			size = 4;
			for(i = 0; i < size; i++){
				for(j = 0; j < size; j++) a[i][j] = 0.;
				if(i > 0) a[i][i - 1] =  1.;
				a[i][i] = 3.;
				if(i < size - 1) a[i][i + 1] = 1.;
				b[i] = 1.;
			}
			break;
	}

	////////// 実行部 //////////

	printf("A :\n"); matprint(a); printf("\n");
	printf("b :\n"); vecprint(b); printf("\n");
	printf("リターンキーを押してください\n"); getchar();

	////////// Jacobi 法 //////////
	printf("----- Jacobi 法 -----\n\n");
	for(i = 0; i < size; i++) x[i] = 0.;
	step = 0;
	e = 1.;	// 誤差の初期値はダミー
	while(e > adm && step < maxrec){
		if(trace){ printf("%d-th step :\n\n", step); vecprint(x); printf("\n"); }
		Jacobi(a, b, x);	// 近似ベクトル x の更新
		matvecmul(a, x, w);	// w = Ax
		vecsub(w, b, c);	// c = w - b = Ax - b
		e = norm(c);		// e = (c の長さ)
		step++;
	}
	printf("%d-th step :\n\n", step); vecprint(x); printf("\n");
	printf("検算 Ax - b :\n\n");
	vecprint(c); printf("\n");
	printf("リターンキーを押してください\n"); getchar();
	
	////////// Gauss-Seidel 法 //////////
	printf("----- Gauss-Seidel 法 -----\n\n");
	for(i = 0; i < size; i++) x[i] = 0.;
	step = 0;
	e = 1.;	// 誤差の初期値はダミー
	while(e > adm && step < maxrec){
		if(trace){ printf("%d-th step :\n\n", step); vecprint(x); printf("\n"); }
		GaussSeidel(a, b, x);	// 近似ベクトル x の更新
		matvecmul(a, x, w);		// w = Ax
		vecsub(w, b, c);		// c = w - b = Ax - b
		e = norm(c);			// e = (c の長さ)
		step++;
	}
	printf("%d-th step :\n\n", step); vecprint(x); printf("\n");
	printf("検算 Ax - b :\n\n");
	vecprint(c); printf("\n");
}

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

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

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

// 行列 a とベクトル v の積をベクトル w に格納する関数
void matvecmul(double a[][N], double v[], double w[])
{
	int i, j;
	double z;
	for(i = 0; i < size; i++){
		z = 0;
		for(j = 0; j < size; j++)
			z += a[i][j] * v[j];
		w[i] = z;
	}
	return;
}

// ベクトル u とベクトル v の差 u - v をベクトル w に格納する関数
void vecsub(double u[], double v[], double w[])
{
	int i;
	for(i = 0; i < size; i++) w[i] = u[i] - v[i];
	return;
}

// ベクトル v の長さを返す関数
double norm(double v[])
{
	int i;
	double z = 0.;
	for(i = 0; i < size; i++) z += v[i] * v[i];
	return sqrt(z);
}