LU分解求线型方程

//用LU分解求线性方程

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

#define N 3
void LUFun();

void main()
{
	LUFun();
	printf("\n");
}


void LUFun()
{
	int i,j,r,k;float m1,m2,y[N],x[N];
	m1=0.0;m2=0.0;
	float l[N][N]={0};
	float u[N][N]={0};
	float a[N][N]={{2,1,1},{1,3,2},{1,2,2}};
	float b[N]={4,6,5};

	for(i=0;i<N;i++)
	{
		u[0][i]=a[0][i];
		l[i][0]=a[i][0]/u[0][0];
		l[i][i]=1;
		/*printf("u[1][%d]=%f\n",i+1,u[0][i]);
		printf("l[%d][1]=%f\n",i+1,l[i][0]);*/
	}

	for(r=1;r<N;r++)
	{
		for(i=r;i<N;i++)
		{
			for(k=0;k<=r-1;k++)
			{
				m1=m1+l[r][k]*u[k][i];
			}

			u[r][i]=a[r][i]-m1;
			m1=0;
			//printf("u[%d][%d]=%f\n",r+1,i+1,u[r][i]);
		}

		for(j=r+1;j<N;j++)
		{
			for(k=0;k<=r-1;k++)
			{
				m2=m2+l[j][k]*u[k][r];
			}
			l[j][r]=(a[j][r]-m2)/u[r][r];
			m2=0;
			//printf("l[%d][%d]=%f\n",j+1,r+1,l[j][r]);
		}
	}
	///////////输出LU//////////////
	printf("L=\n");
	for(i=0;i<N;i++)
	{
		for(j=0;j<N;j++)
		{
			printf("%-10f",l[i][j]);
		}
		printf("\n");
	}
	printf("\nU=\n");
	for(i=0;i<N;i++)
	{
		for(j=0;j<N;j++)
		{
			printf("%-10f",u[i][j]);
		}
		printf("\n");
	}

	/****解方程************/
	m1=0;m2=0;
	y[0]=b[0];
	printf("y1=%f\t",y[0]);
	for(i=1;i<N;i++)
	{
		for(k=0;k<=i-1;k++)
		{m1=m1+l[i][k]*y[k];}
		y[i]=b[i]-m1;m1=0;
		printf("y%d=%f\t",i+1,y[i]);
	}

	x[N-1]=y[N-1]/u[N-1][N-1];
	printf("\nx%d=%f\t",N,x[N-1]);
	for(i=N-2;i>=0;i--)
	{
		for(k=i+1;k<N;k++)
		{m2=m2+u[i][k]*x[k];}
		x[i]=(y[i]-m2)/u[i][i];m2=0;
		printf("x%d=%f\t",i+1,x[i]);
	}
}