//////////////////////////////////////////////////////////
//
// LU.c  数値解析(C) 2008.12.24 の課題
// LU 分解法
// ３次元限定のサンプルプログラム
//
//////////////////////////////////////////////////////////
//
// 課題：
// ・大域変数 SIZE をサイズとする行列が扱えるように、
// 　LU_decomp, LU_step2, LU_step3 を書き換えよ
// ・main() 内で SIZE の値を 4, 8, 16 と取り替えて実行せよ。
//
// 発展課題：
// ・行列やベクトルを取り替えて実行してみよ。
//
//////////////////////////////////////////////////////////

#include"sqmatrix.h"
// sqmatrix.h を同じディレクトリに download して用いよ。
// ( sqmatrix.h は書き換えなくて良い。)

// 正方行列  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; 

	// 課題の行列・ベクトル
	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");
}

/*　実行例 ( SIZE = 3 のとき )

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
*/