LU算法 源程序

#include<math.h>
#include<stdio.h>
#define N 3
static double aa[N][N]={{3,1,6},{2,1,3},{1,1,1}};
static double bb[N]={2,7,4};
int LU(double a[][N+1])//返回：-1存在a[k,k]=0,1正常结束{
int i,r,k;
double add1,add2;
for(i=1;i<=N;i++)//{//第一行U[1][i]
a[1][i]=a[1][i];
//printf("%12lf",a[1][i]);}
for(i=2;i<=N;i++)//第一列L[i][1]
{
if(a[1][1]==0)
return -1;
a[i][1]=a[i][1]/a[1][1];
//printf("%12lf",a[i][1]);
}
for(r=2;r<=N;r++)//第r行，第r列
{
for(i=r;i<=N;i++)//行U[][]
{
add1=0;
for(k=1;k<=r-1;k++)
{
add1=add1+a[r][k]*a[k][i];
}
a[r][i]=a[r][i]-add1;//printf("%12lf",a[r][i]);
}
for(i=r+1;i<=N;i++)//列L[][]
{
add2=0;
for(k=1;k<=r-1;k++)
{
add2=add2+a[i][k]*a[k][r];
}
if(a[r][r]==0)
return -1;
a[i][r]=(a[i][r]-add2)/a[r][r];//printf("%12lf",a[i][r]);
}
}
return 1;
}//对N阶矩阵A进行LU分解，结果L、U存放在A的相应位置
void solve(double a[][N+1],double b[N+1])//利用分解的LU求x,,,LUx=b
{
int i,k;
double add3,add4;
/*for(i=1;i<=N;i++)
printf("%12lf",b[i]);*/
b[1]=b[1];//printf("%12lf",b[1]);
for(i=2;i<=N;i++)//第一次回代求y[i];;;Ly=b
{
add3=0;
for(k=1;k<i;k++)
add3=add3+a[i][k]*b[k];
b[i]=b[i]-add3;//printf("%12lf",b[i]);
}
/*if(a[N][N]==0)
return -1;*/
b[N]=b[N]/a[N][N];//printf("%12lf",b[N]);
for(i=N-1;i>0;i--)//第二次回代求x[i];;;Ux=y
{
add4=0;
/*if(a[i][i]==0)
return -1;*/
for(k=i+1;k<=N;k++)
add4=add4+a[i][k]*b[k];
b[i]=(b[i]-add4)/a[i][i];//printf("%12lf",b[i]);
}
}//回代求解，L与U元素均在A矩阵相应位置；b存放常数列，返回时存放方程组的解void main()
{
int d,i,j;
double a[N+1][N+1],b[N+1];
for(i=1;i<=N;i++)//初始化
{
for(j=1;j<=N;j++)
a[i][j]=aa[i-1][j-1];
b[i]=bb[i-1];
}
d=LU(a);
if(d==1)
{
printf("矩阵L如下\n");
for(i=1;i<=N;i++)
{
for(j=1;j<=i;j++)
if(i==j)
printf("%12d ",1);
else
printf("%12lf",a[i][j]);
printf("\n");
}
printf("\n矩阵U如下\n");
for(i=1;i<=N;i++)
{
for(j=1;j<=N;j++)
if(i<=j)
printf("%12lf",a[i][j]);
else
printf("%12s","");
printf("\n");
}
solve(a,b);
for(i=1;i<=N;i++)
printf(" x[%d]=%lf ",i,b[i]);
printf("\n");
}
}。