LU-decomposition.rar_decomposition

#include<stdio.h> #define N 100 /*输入函数*/ void input(float a[][N],int n) { int i,j; for(i=0;i<n;i++) for(j=0;j<n+1;j++) scanf("%f",&a[i][j]); } /*拆分函数*/ void split(float A[][N],float a[][N],float b[][1],int n) { int i,j; for(i=0;i<n;i++) for(j=0;j<n;j++) a[i][j]=A[i][j]; for(i=0;i<n;i++) b[i][0]=A[i][n]; } /*输出函数*/ void output(float x[][1],int n) { int i; for(i=0;i<n;i++) printf("x(%d)=%f\n",i+1,x[i][0]); } /*消元函数*/ void eliminate(float a[][N],float b[][1],int k,int n) { int i,j; float m; for(i=k+1;i<n;i++) { m=a[i][k]/a[k][k]; b[i][0]=b[i][0]-m*b[k][0]; for(j=k;j<n;j++) a[i][j]=a[i][j]-m*a[k][j]; } } /*回代函数*/ void back(float a[][N],float b[][1],int n) { int i,j; float s; b[n-1][0]=b[n-1][0]/a[n-1][n-1]; for(i=n-2;i>=0;i--) { s=0; for(j=i+1;j<n;j++) s=s+a[i][j]*b[j][0]; b[i][0]=(b[i][0]-s)/a[i][i]; } } /*初始化函数*/ void initialize(float a[N][N],int n) { int i,j; for(i=0;i<n;i++) for(j=0;j<n;j++) a[i][j]=0; } /*矩阵分解函数*/ void disassemble(float a[][N],float b[][1],float l[][N],float u[][N],int n) { int i,j,k; float s,m,p; for(i=0;i<n;i++) u[0][i]=a[0][i]; for(i=0;i<n;i++) l[i][0]=a[i][0]/u[0][0]; for(i=1;i<n;i++) { for(j=i;j<n;j++) { s=0; m=0; for(k=0;k<i;k++) s+=l[i][k]*u[k][j]; u[i][j]=a[i][j]-s; for(k=0;k<i;k++) m+=l[j][k]*u[k][i]; l[j][i]=(a[j][i]-m)/u[i][i]; } p=0; for(k=0;k<i;k++) p+=l[i][k]*b[k][0]; b[i][0]=b[i][0]-p; } } /*显示函数*/ void show(float a[][N],int n) { int i,j; for(i=0;i<n;i++) { for(j=0;j<n;j++) printf("%f ",a[i][j]); printf("\n"); } } /*主函数*/ void main() { float A[N][N],a[N][N],b[N][1],u[N][N],l[N][N]; int n; printf("请输入方程未知数的个数：\n"); scanf("%d",&n); initialize(u,n); /*初始化u矩阵*/ initialize(l,n); /*初始化l矩阵*/ printf("请输入所求线性代数方程组的增广矩阵：\n"); input(A,n); /*输入矩阵*/ split(A,a,b,n); /*将增广矩阵拆分为系数矩阵与列矩阵*/ disassemble(a,b,l,u,n); /*拆分系数阵*/ printf("L=\n"); show(l,n); /*显示l矩阵*/ printf("U=\n"); show(u,n); /*显示u矩阵*/ back(u,b,n); /*回代*/ printf("求解方程为：\n"); output(b,n); /*输出结果*/ }