/* lu.c  数値解析Ｃ 2003.10.29 の課題
LU 分解法
gcc lu.c ; a.out
*/

#include<stdio.h>
#include<stdlib.h>
#define RANDOMIZE()	srand(time(NULL))

/* 扱うことのできるベクトル・行列の最大サイズ設定 */

#define SIZE 128

/* 縦ベクトルの型定義 
   row : サイズ、
   ent : 成分 
を表す */

typedef struct {
	int row;
	double ent[SIZE];
} VECTOR;

/* 行列の型定義 
   row : 行数、
   col : 列数、
   ent : 成分 
を表す */

typedef struct {
	int row,col;
	double ent[SIZE][SIZE]; 
} MATRIX;

/* ベクトルを出力する手続き */

void vectorwrite(VECTOR a)
{
	int i;
	int m=a.row;

	for(i=1;i<=m;i++) printf("%12.7lf\n",a.ent[i]);
}

/* 行列を出力する手続き */

void matrixwrite(MATRIX a)
{
	int i,j;
	int m=a.row,n=a.col;
	
	for(i=1;i<=m;i++){
		for(j=1;j<=n;j++)
			printf("%12.7lf",a.ent[i][j]);
		printf("\n");
	}
}
 
/* 行列の積  C = A B  をを返す関数
   A の列数 != B の行数  のときはエラーメッセージを出力する */
   
MATRIX matrixmult(MATRIX a, MATRIX b)
{
	int i,j,k;
	int l=a.row,m=a.col,n=b.col;
	MATRIX c;
	double x;

	if(a.col!=b.row) printf("Invalid matrix sizes.\n");
	c.row=l;
	c.col=n;
	for(i=1;i<=l;i++){
		for(j=1;j<=n;j++){
			x=0.0;
			for(k=1;k<=m;k++)
				x+=a.ent[i][k]*b.ent[k][j];
			c.ent[i][j]=x;
		}
	}
	return c;
}

/* 行列とベクトルの積  c = A b  を返す関数
   A の列数 != b の項数  のときはエラーメッセージを出力する */
   
VECTOR matvecmult(MATRIX a, VECTOR b)
{
	int i,k;
	int l=a.row,m=a.col;
	VECTOR c;
	double x;

	if(a.col!=b.row) printf("Invalid matrix sizes.\n");
	c.row=m;
	for(i=1;i<=l;i++){
		x=0.0;
		for(k=1;k<=m;k++)
			x+=a.ent[i][k]*b.ent[k];
		c.ent[i]=x;
	}
	return c;
}

/* 正方行列  A  の LU 分解を計算する手続き
   ふたつの行列  L, U  はポインタで返す
   A  が正方行列でないときはエラーメッセージを出力する
   pivot 選択は使っていない */

void lu_decomp(MATRIX a, MATRIX *l, MATRIX *u)
{
	int i,j,k,m,n=a.row;
	double x;

	if(a.row!=a.col) printf("Invalid matrix sizes.\n");
	l->row=n;
	l->col=n;
	u->row=n;
	u->col=n;
	for(i=1;i<=n;i++){
		for(j=1;j<=n;j++){
			l->ent[i][j]=0.0;
			u->ent[i][j]=0.0;
		}
	}
	for(i=1;i<=n;i++)
		u->ent[i][i]=1.0;

	for(k=1;k<=n;k++){
		for(j=1;j<=k;j++){
			x=0.0;
			for(m=1;m<j;m++) 
				x+=l->ent[k][m]*u->ent[m][j];
			l->ent[k][j]=a.ent[k][j]-x;
		}
		for(j=k+1;j<=n;j++){
			x=0.0;
			for(m=1;m<j;m++) 
				x+=l->ent[k][m]*u->ent[m][j];
			u->ent[k][j]=(a.ent[k][j]-x)/l->ent[k][k];
		}
	}
}

/* LU 分解法の前進代入 :  L y = b  の解  y  を求める */

VECTOR step2(MATRIX l, VECTOR b)
{
	int i,j;
	int n=l.row;
	VECTOR y;
	double w;

	y.row=n;
	for(i=1;i<=n;i++){
		w=0.0;
		for(j=1;j<i;j++)
			w+=l.ent[i][j]*y.ent[j];
		y.ent[i]=(b.ent[i]-w)/l.ent[i][i];
	}
	return y;
}

/* LU 分解法の後退代入 :  U x = y  の解  x  を求める */

VECTOR step3(MATRIX u, VECTOR y)
{
	int i,j;
	int n=u.row;
	VECTOR x;
	double w;

	x.row=n;
	for(i=n;i>=1;i--){
		w=0.0;
		for(j=i+1;j<=n;j++)
			w+=u.ent[i][j]*x.ent[j];
		x.ent[i]=(y.ent[i]-w)/u.ent[i][i];
	}
	return x;
}

/* ランダムな (m,n) 行列を返す関数 */

MATRIX randommatrix(int m, int n)
{
	int i,j;
	MATRIX a;
	double x;

	a.row=m;
	a.col=n;
	x=rand();
	for(i=1;i<=m;i++)
		for(j=1;j<=n;j++)
			a.ent[i][j]=rand()/x;
	return a;
}

main()
{
	MATRIX a,l,u,c;
	VECTOR b,x,y,z;
	int i,j,n,k;

	RANDOMIZE();

	for(k=1;k<=9;k++){
		printf("----------\n");
		printf("\n");
		n=5;

/* 課題の行列  A  の定義 */

		a.row=n;
		a.col=n;
		for(i=1;i<=n;i++) for(j=1;j<=n;j++) a.ent[i][j]=1.0;
		for(i=1;i<=n;i++) a.ent[i][i]=k/10.0;
		a.ent[1][n]=2.0;
		a.ent[n][1]=2.0;

/* 乱数行列で実験するときはこちらを使う */
/*
		a=randommatrix(n,n);
*/

/* 課題のベクトル  b  の定義 */

		b.row=n;
		for(i=1;i<=n;i++) b.ent[i]=i;
		
		lu_decomp(a,&l,&u);
		printf("A :\n"); matrixwrite(a); printf("\n");
		printf("L :\n"); matrixwrite(l); printf("\n");
		printf("U :\n"); matrixwrite(u); printf("\n");
		c=matrixmult(l,u);
		printf("検算 LU :\n"); matrixwrite(c); printf("\n");
		y=step2(l,b);
		x=step3(u,y);
		z=matvecmult(a,x);
		printf("b :\n"); vectorwrite(b); printf("\n");
		
/* y  も出力するときはコメントアウトをはずす	
		printf("y :\n"); vectorwrite(y); printf("\n"); 
*/
		printf("x :\n"); vectorwrite(x); printf("\n");
		printf("検算 Ax :\n"); vectorwrite(z); printf("\n");
	}
}

/*　実行例

----------

A :
   0.8991140   1.0011075   0.9288434   0.3419395   0.6698969
   1.0094137   0.9597840   0.2149581   0.8481346   0.8137330
   0.9329619   0.2217069   0.2892296   0.7064096   0.1779608
   0.9617222   0.2282481   0.7629958   0.4179414   0.3855126
   1.0266837   0.4069703   0.3443275   1.0260608   0.6266353

L :
   0.8991140   0.0000000   0.0000000   0.0000000   0.0000000
   1.0094137  -0.1641353   0.0000000   0.0000000   0.0000000
   0.9329619  -0.8170881   3.4464799   0.0000000   0.0000000
   0.9617222  -0.8425697   4.0190540  -0.0459436   0.0000000
   1.0266837  -0.7361782   2.9966804   0.2571244   2.0559710

U :
   1.0000000   1.1134378   1.0330652   0.3803072   0.7450633
   0.0000000   1.0000000   5.0435944  -2.8284422  -0.3756401
   0.0000000   0.0000000   1.0000000  -0.5685478  -0.2391094
   0.0000000   0.0000000   0.0000000   1.0000000  -6.8227029
   0.0000000   0.0000000   0.0000000   0.0000000   1.0000000

検算 LU :
   0.8991140   1.0011075   0.9288434   0.3419395   0.6698969
   1.0094137   0.9597840   0.2149581   0.8481346   0.8137330
   0.9329619   0.2217069   0.2892296   0.7064096   0.1779608
   0.9617222   0.2282481   0.7629958   0.4179414   0.3855126
   1.0266837   0.4069703   0.3443275   1.0260608   0.6266353

b :
   1.0000000
   2.0000000
   3.0000000
   4.0000000
   5.0000000

x :
   1.0300851
  -5.5297696
   1.9024030
   2.7517483
   4.3316580

検算 Ax :
   1.0000000
   2.0000000
   3.0000000
   4.0000000
   5.0000000

*/