- #include <stdio.h>
- #include <stdlib.h>
- //构造动态二维数组
- float** create_Array(int n, int m)
- {
- float** a = (float**) malloc(sizeof(float*)*n);
- for (int k=0; k<n; ++k)
- a[k] = (float*)malloc(sizeof(float)*m);
- return a;
- }
- //打印二维数组
- void print_Array(float** a,int n,int m)
- {
- for(int i=0;i<n;i++)
- {
- for(int j=0;j<m;j++)
- printf("%0.5f\\t\\t",a[i][j]);
- printf("\\n");
- }
- }
- //LU分解法回代过程,并求出方程组的解
- void get_X(float** &L,float** &U,float* &b,int n)
- {
- int i,j;
- float* y=(float*)malloc(sizeof(float)*n);
- float* x=(float*)malloc(sizeof(float)*n);
- //计算Ly=b中的y
- for(i=0;i<n;i++)
- {
- y[i]=b[i];
- for(j=0;j<=i-1;j++)
- y[i]-=(L[i][j]*y[j]);
- y[i]/=L[i][i];
- }
- //计算Ux=y中的x
- for(i=n-1;i>=0;i--)
- {
- x[i]=y[i];
- for(j=i+1;j<n;j++)
- x[i]-=(U[i][j]*x[j]);
- x[i]/=U[i][i];
- }
- //输出x
- for(i=0;i<n;i++)
- printf("x[%d]=%f\\n",i+1,x[i]);
- free(y);
- free(x);
- }
- //释放动态二维数组的内存空间
- void free_Array(float **a,int n)
- {
- for(int m=0;m <n;m++)
- delete[] a[m];
- delete[] a;
- }
- void main()
- {
- int n,m;
- int i,j,k;
- float sum;
- //建立并输入增广矩阵
- printf("Input row and rank of matrix:\\n");
- scanf("%d%d",&n,&m);
- float** a=create_Array(n,m);
- printf("Input the matrix:\\n");
- for(i=0;i<n;i++)
- for(j=0;j<m;j++)
- {
- scanf("%f",&a[i][j]);
- }
- float *b=(float*)malloc(sizeof(float)*n);
- for(i=0;i<n;i++)
- b[i]=a[i][m-1];
- //打印增广矩阵
- printf("a:\\n");
- print_Array(a,n,m);
- float** L=create_Array(n,n);
- float** U=create_Array(n,n);
- //初始化L矩阵
- for(i=0;i<n;i++)
- for(j=0;j<n;j++)
- {
- if(i==j)
- L[i][j]=1;
- else
- L[i][j]=0;
- }
- //初始化U矩阵
- for(i=0; i<n; i++)
- {
- U[0][i] = (float)(a[0][i]/L[0][0]);
- }
- for(i=1;i<n;i++)
- for(j=0;j<n;j++)
- U[i][j]=0;
- //计算出L和U矩阵
- for(i=0; i<n-1; i++)
- {
- for(j=i+1; j<n; j++)
- {
- for(k=0,sum=0; k<n; k++)
- {
- if(k != i)
- sum += L[j][k]*U[k][i];
- }
- L[j][i] = (float)((a[j][i]-sum)/U[i][i]);
- }
- for(j=i+1; j<n; j++)
- {
- for(k=0,sum=0; k<n; k++)
- {
- if(k != i+1)
- sum += L[i+1][k]*U[k][j];
- }
- U[i+1][j] = (float)((a[i+1][j]-sum));
- }
- }
- //打印L矩阵
- printf("L:\\n");
- print_Array(L,n,n);
- //打印U矩阵
- printf("U:\\n");
- print_Array(U,n,n);
- //求解
- get_X(L,U,b,n);
- //释放内存
- free_Array(a,n);
- free_Array(U,n);
- free_Array(L,n);
- free(b);
- }
- //该片段来自于http://www.codesnippet.cn/detail/080520133200.html
来源: http://www.codesnippet.cn/detail/080520133200.html