/* JGS.c
連立一次方程式の反復解法: Jacobi 法, Gauss-Seidel 法
K.Shiota
2006.1-2008.1
*/

#include"sqmatrix.h"
// sqmatrix.h を同じフォルダにダウンロードせよ

VECTOR JacobiStep(MATRIX a, VECTOR b, VECTOR x)
{
	VECTOR w;
	int i,j;
	double z;
	for(i=0;i<SIZE;i++){
		z=b.ent[i];
		for(j=0;j<SIZE;j++)
			if(j!=i) z-=a.ent[i][j]*x.ent[j];
		w.ent[i]=z/a.ent[i][i];
	}
	return w;
}

VECTOR GSStep(MATRIX a, VECTOR b, VECTOR x)
{
	VECTOR w;
	int i,j;
	double z;
	for(i=0;i<SIZE;i++){
		z=b.ent[i];
		for(j=0;j<i;j++)      z-=a.ent[i][j]*w.ent[j];
		for(j=i+1;j<SIZE;j++) z-=a.ent[i][j]*x.ent[j];
		w.ent[i]=z/a.ent[i][i];
	}
	return w;
}

double err(MATRIX a, VECTOR b, VECTOR x)
{
	return vecnorm(vecsub(matvecmul(a,x),b));
}

main()
{
	MATRIX a;
	VECTOR b,x;
	int i,maxrec=128;
	double adm=1e-7;

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

/* テストデータ */

/*
	// 収束する例
	SIZE=2;
	a.ent[0][0]=3; a.ent[0][1]=-2; a.ent[1][0]=1; a.ent[1][1]=3;
	b.ent[0]=1; b.ent[1]=4;
*/

/*
	// 収束しない例
	SIZE=2;
	a.ent[0][0]=1; a.ent[0][1]=3; a.ent[1][0]=3; a.ent[1][1]=-2;
	b.ent[0]=4; b.ent[1]=1;
*/

/*
	// 式・未知数を入れ替えることにより収束の速さが変わる例
	SIZE=2;
	a.ent[0][0]=5; a.ent[0][1]=4;
	a.ent[1][0]=-2; a.ent[1][1]=3;
	b.ent[0]=9; b.ent[1]=1;
*/

/*
	SIZE=2;
	a.ent[0][0]=3; a.ent[0][1]=-2;
	a.ent[1][0]=4; a.ent[1][1]=5;
	b.ent[0]=1; b.ent[1]=9;
*/

/*
	// 式・未知数を入れ替えることにより収束の速さが変わる例
	SIZE=3;
	a.ent[0][0]=3; a.ent[0][1]=2; a.ent[0][2]=2;
	a.ent[1][0]=2; a.ent[1][1]=5; a.ent[1][2]=-4;
	a.ent[2][0]=1; a.ent[2][1]=1; a.ent[2][2]=7;
	b.ent[0]=7; b.ent[1]=3; b.ent[2]=9;
*/
/*
	SIZE=3;
	a.ent[0][0]=7; a.ent[0][1]=1; a.ent[0][2]=1;
	a.ent[1][0]=-4; a.ent[1][1]=5; a.ent[1][2]=2;
	a.ent[2][0]=2; a.ent[2][1]=2; a.ent[2][2]=3;
	b.ent[0]=9; b.ent[1]=3; b.ent[2]=7;
*/

	// ランダムな整数係数の場合
	SIZE=8;
	a=randomintmat(6);
	for(i=0;i<SIZE;i++)
		a.ent[i][i]=RANDOM(6*SIZE)+1;
	b=randomintvec(6);

/*
	// ランダムな係数の場合
	SIZE=64;
	a=randommat(5);
	for(i=0;i<SIZE;i++)
		a.ent[i][i]*=SIZE;
	b=randomvec(5);
*/

/* 実行部 */

	printf("A :¥n"); matprint(a); printf("¥n");
	printf("b :¥n"); vecprint(b); printf("¥n");
	printf("Enter Key を押してください。¥n");
	getchar();

	/* Jacobi 法 */
	printf("----- Jacobi 法 -----¥n¥n");
	x=zerovec();
	i=0;
	// printf("%d-th step :¥n¥n",i); vecprint(x); printf("¥n");
	while(err(a,b,x)>adm && i<maxrec){
		x=JacobiStep(a,b,x);
		i++;
		// printf("%d-th step :¥n¥n",i); vecprint(x); printf("¥n");
	}
	printf("%d-th step :¥n¥n",i); vecprint(x); printf("¥n");
	printf("検算 Ax-b :¥n¥n");
	vecprint(vecsub(matvecmul(a,x),b)); printf("¥n");
	printf("Enter Key を押してください。¥n");
	getchar();
	
	/* GS 法 */
	printf("----- Gauss-Seidel 法 -----¥n¥n");
	x=zerovec();
	i=0;
 	printf("%d-th step :¥n¥n",i); vecprint(x); printf("¥n");
	while(err(a,b,x)>adm && i<maxrec){
		x=GSStep(a,b,x);
		i++;
		printf("%d-th step :¥n¥n",i); vecprint(x); printf("¥n");
	}
	printf("%d-th step :¥n¥n",i); vecprint(x); printf("¥n");
	printf("検算 Ax-b :¥n¥n");
	vecprint(vecsub(matvecmul(a,x),b)); printf("¥n");
}