列主元lu分解matlab程序,矩阵的列选主元LU分解程序.请大家帮忙精简一下，谢谢！...

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]);

*/

}