C语言程序如下:
#include
#define MAX 20
main()
{int i,j,n;
int A[MAX][MAX];
void devide(int A[][MAX],int n);//函数声明
printf("请输入方程组系数数组的维数n=");
scanf("%d",&n);
putchar('\n');putchar('\n');
putchar('\n');putchar('\n');
printf("请输入方程组系数(按行):");
putchar('\n');
for(i=1;i<=n;i++)
{for(j=1;j<=n;j++)
scanf("%d",&A[i][j]);
}
printf("方程组系数矩阵A如下:");
putchar('\n');
for(i=1;i<=n;i++)
{for(j=1;j<=n;j++)
printf("%d ",A[i][j]);
putchar('\n');
putchar('\n');putchar('\n');putchar('\n');
}
devide(A,n);//分解函数调用
return 0;
}
void devide(int A[][MAX],int n)//定义列选主元LU分解函数
{int i,j;
int k=1;//控制在进行第几框架的分解
float M=0;
int K;
float max;
int line;
int cow;
float temp;
float U[MAX][MAX];
float L[MAX][MAX];
int P[MAX][MAX];
float y[MAX][1];
float x[MAX][1];
int B[MAX][1];
printf("请输入方程组等号右侧常数(按行):");
putchar('\n');
putchar('\n');
putchar('\n');
for(i=1;i<=n;i++)
{scanf("%d",&B[i][1]);
}
for(i=1;i<=n;i++)//为L,U,P赋初值为0
{for(j=1;j<=n;j++)
{U[i][j]=0;
L[i][j]=0;
if(i==j)
P[i][j]=1;
else
P[i][j]=0;
}
}
if(k==1)//如果当前正在进行第一框架的分解
{max=A[1][1];
line=1;
for(i=2;i<=n;i++)//当k=1时,选主元的情况
{if(A[i][1]>max)
{max=A[i][1];
line=i;}
}
for(cow=1;cow<=n;cow++)//交换两行
{temp=A[1][cow];
A[1][cow]=A[line][cow];
A[line][cow]=temp;
}
temp=B[1][1];
B[1][1]=B[line][1];
B[line][1]=temp;
for(i=1;i<=n;i++)
{ U[1][i]=A[1][i];//A的第一行为U的第一行
}//L的第一列为A的第一列除以U对角线上的元素
for(i=1;i<=n;i++)
{L[i][1]=A[i][1]/U[1][1];}
k=2;
}
while((k>=2)&&(k<=n))//如果当前正在进行第二框架以上的分解
{//先选取最大的列元
for(j=1;j<=k-1;j++)
{M+=L[k][j]*U[j][k];}
U[k][k]=A[k][k]-M;
max=U[k][k];
M=0;
line=k;
for(j=k;j<=n;j++)
{for(K=1;K<=k-1;K++)
{M+=L[j][K]*U[K][k];}
U[j][k]=A[j][k]-M;
M=0;
if(U[j][k]>max)
{max=U[j][k];
line=j;}
U[j][k]=0;
}
for(cow=1;cow<=n;cow++)//交换两行
{temp=A[k][cow];
A[k][cow]=A[line][cow];
A[line][cow]=temp;
temp=L[k][cow];
L[k][cow]=L[line][cow];
L[line][cow]=temp;
}
temp=B[1][1];
B[1][1]=B[line][1];
B[line][1]=temp;
for(j=k;j<=n;j++)
{for(K=1;K<=k-1;K++)
{M+=L[k][K]*U[K][j];}
U[k][j]=A[k][j]-M;
M=0;}
for(j=k;j<=n;j++)
{for(K=1;K<=k-1;K++)
{ M+=L[j][K]*U[K][k];}
L[j][k]=(A[j][k]-M)/U[k][k];
M=0;}
k++;
}
printf("分解后U为:");
putchar('\n');
putchar('\n');
for(i=1;i<=n;i++)//显示输出分解后的L,U,P
{for(j=1;j<=n;j++)
printf("%f ",U[i][j]);
putchar('\n');putchar('\n');putchar('\n');}
putchar('\n');
putchar('\n');
putchar('\n');
putchar('\n');
printf("分解后L为:");
putchar('\n');
putchar('\n');
for(i=1;i<=n;i++)
{for(j=1;j<=n;j++)
printf("%f ",L[i][j]);
putchar('\n');putchar('\n');putchar('\n');}
putchar('\n');
putchar('\n');
putchar('\n');
putchar('\n');
M=0;
//求方程组的解:LUx=b,Ux=y
for(i=1;i<=n;i++)
{if(i==1)
y[i][1]=B[1][1];
else
{
for(j=1;j<=i-1;j++)
{M+=L[i][j]*y[j][1];}
y[i][1]=B[i][1]-M;
M=0;}
}
for(i=n;i>=1;i--)//求xn
{if(i==n)
x[i][1]=y[i][1]/U[i][i];
else
{
for(j=n;j>=i+1;j--)
{M+=U[i][j]*x[j][1];}
x[i][1]=(y[i][1]-M)/U[i][i];
M=0;}
}
printf("方程组的解为:");
putchar('\n');
for(i=1;i<=n;i++)
{printf("x=%f ",x[i][1]);
putchar('\n');
}
/*printf("%f ",y[1][1]);
printf("%f ",y[2][1]);
printf("%f ",y[3][1]);
*/
}