/* LU.c  数値解析Ｃ 2007.12.19 の課題
LU 分解法
３次元限定version
gcc LU.c -lm
*/

#include"sqmatrix.h"
// sqmatrix.h を同じディレクトリに download せよ

/* 正方行列  A  の LU 分解を計算する手続き
   ふたつの行列  L, U  はポインタで返す
   pivot 選択は使っていない */
/***** ３次元限定なのでここを書き換えよ *****/
void LU_decomp(MATRIX A, MATRIX *L, MATRIX *U)
{
	*L=zeromat(); // 初期値はゼロ行列
	*U=unitmat(); // 初期値は単位行列

	L->ent[0][0] = A.ent[0][0];
	U->ent[0][1] = A.ent[0][1] / L->ent[0][0];
	U->ent[0][2] = A.ent[0][2] / L->ent[0][0];

	L->ent[1][0] = A.ent[1][0];
	L->ent[1][1] = A.ent[1][1] - L->ent[1][0] * U->ent[0][1];
	U->ent[1][2] = (A.ent[1][2] - L->ent[1][0] * U->ent[0][2]) / L->ent[1][1];

	L->ent[2][0] = A.ent[2][0];
	L->ent[2][1] = A.ent[2][1] - L->ent[2][0] * U->ent[0][1];
	L->ent[2][2] = A.ent[2][2] - L->ent[2][0] * U->ent[0][2]
                                   - L->ent[2][1] * U->ent[1][2];
}

/* LU 分解法の前進代入 :  L y = b  の解  y  を求める */
/***** ３次元限定なのでここも書き換えよ *****/
VECTOR LU_step2(MATRIX L, VECTOR b)
{
	VECTOR y;
	y.ent[0] = b.ent[0] / L.ent[0][0];
	y.ent[1] = (b.ent[1] - L.ent[1][0] * y.ent[0]) / L.ent[1][1];
	y.ent[2] = (b.ent[2] - L.ent[2][0] * y.ent[0] 
                             - L.ent[2][1] * y.ent[1]) / L.ent[2][2];
	return y;
}

/* LU 分解法の後退代入 :  U x = y  の解  x  を求める */
/***** ３次元限定なのでここも書き換えよ *****/
VECTOR LU_step3(MATRIX U, VECTOR y)
{
	VECTOR x;
	x.ent[2] = y.ent[2];
	x.ent[1] = y.ent[1] - U.ent[1][2] * x.ent[2];
	x.ent[0] = y.ent[0] - U.ent[0][1] * x.ent[1]
                            - U.ent[0][2] * x.ent[2];
	return x;
}

main()
{
	MATRIX A,L,U,C;
	VECTOR b,x,y,z;
	int i,j,n,k;

	RANDOMIZE();
	SIZE = 3; // 行列のサイズを設定

	/***** 課題 (3) の行列・ベクトル *****/
	for(i=0; i<SIZE; i++)
		for(j=0; j<SIZE; j++)
			A.ent[i][j] = 1.0/(i+j+1);
	for(i=0; i<SIZE; i++)
		b.ent[i] = i+1;

	/* 乱数行列・乱数ベクトルで実験するとき */
//	A = randommat(10.0);
//	b = randomvec(10.0);

	LU_decomp(A,&L,&U);
	printf("A :\n"); matprint(A); printf("\n");
	printf("L :\n"); matprint(L); printf("\n");
	printf("U :\n"); matprint(U); printf("\n");
	C = matmul(L,U);
	printf("検算 LU :\n"); matprint(C); printf("\n");
	y = LU_step2(L,b);
	x = LU_step3(U,y);
	z = matvecmul(A,x);
	printf("b :\n"); vecprint(b); printf("\n");		
	/* y を出力するときはコメントアウトをはずす */
	//	printf("y :\n"); vecprint(y); printf("\n"); 
	printf("x :\n"); vecprint(x); printf("\n");
	printf("検算 Ax :\n"); vecprint(z); printf("\n");
}

/*　実行例

A :
  -1.6994587  -3.7512074   0.2549254
  -0.1404254   1.3622057   3.2961768
  -2.2395422   3.1518807   3.0007095

L :
  -1.6994587   0.0000000   0.0000000
  -0.1404254   1.6721660   0.0000000
  -2.2395422   8.0952118 -13.1905515

U :
   1.0000000   2.2072954  -0.1500039
   0.0000000   1.0000000   1.9586048
   0.0000000   0.0000000   1.0000000

検算 LU :
  -1.6994587  -3.7512074   0.2549254
  -0.1404254   1.3622057   3.2961768
  -2.2395422   3.1518807   3.0007095

b :
   0.9990628
  -1.8782601
  -1.7182739

x :
  -0.1895321
  -0.2137354
  -0.4895744

検算 Ax :
   0.9990628
  -1.8782601
  -1.7182739
*/