#include <iostream>
#include <math.h>
using namespace std;
void daolite(double a[3][3])
{ //由a[][]分解为L[][](单位下三角矩阵),U[][](上三角矩阵)
double u[3][3],double y[3],double x[3];
double l[3][3];
double b[3]={1.0,3.0,4.0};
int i,r,j,k;
for(r=0;r<3;r++)
for(j=0;j<3;j++)
{
u[r][j]=0;
}
for(r=0;r<3;r++)
for(j=0;j<3;j++)
{
if(r==j) l[r][j]=1;
else l[r][j]=0;
}
for( r=0;r<3;r++)
{
for( j=r;j<3;j++)
{
double s=0.0;
for( k=0;k<r;k++)
s+=l[r][k]*u[k][j];//计算：∑l[r][k]*u[k][j]
u[r][j]=a[r][j]-s;//计算：u[r][j]=a[r][j]-∑l[r][k]*u[k][j]
}
for( i=r+1;i<3;i++)
{
double s=0.0;
for( k=0;k<r;k++)
s+=l[i][k]*u[k][r];//计算：∑l[i][k]*u[k][r]
l[i][r]=(a[i][r]-s)/u[r][r];//计算：l[i][r]=(a[i][r]-∑l[i][k]*u[k][r])/u[r][r]
}
}
cout<<"其中u矩阵为："<<endl;//输出u矩阵
for( r=0;r<3;r++)
{
for( j=0;j<3;j++)
{
cout<<u[r][j]<<"\t";
}
cout<<endl;
}
cout<<"其中l矩阵为："<<endl;//输出l矩阵
for( i=0;i<3;i++)
{
for( r=0;r<3;r++)
{
cout<<l[i][r]<<"\t";
}
cout<<endl;
}
for( i=0;i<3;i++)
{double s=0.0;
for(k=0;k<i;k++)
{
s=s+l[i][k]*y[k];
}
y[i]=b[i]-s;
}
for(i=2;i>=0;i--)
{ double sum=0.0;
for(k=i+1;k<2;k++)
{
sum+=u[i][k]*x[k];
}
x[i]=(y[i]-sum)/u[i][i];
}
cout<<"方程式组的解为：x=(";//输出x[]
for (i=0;i<3;i++)
{
cout<<x[i];
if(i!=2)cout<<",";
}
cout<<")"<<endl;
}
int main()
{
double a[3][3];
int i,j;
cout<<"请输入矩阵:"<<endl;
for(i=0;i<3;i++)
for(j=0;j<3;j++)
cin>>a[i][j];
cout<<"你输入的矩阵为:"<<endl;
for( i=0;i<3;i++)
{
for(j=0;j<3;j++)
{
cout<<a[i][j]<<"\t";
}
cout<<endl;
}
daolite(a);
return 0;
}
