LU矩阵分解法解线性方程组（C++实现）

LU矩阵分解法解线性方程组（C++实现）

#include<cstdio>
#include<iostream>
using namespace std;
void fenjie(double a[],double b[],int n){
	double aa[n][n];
	for(int i=0;i<n;i++){
		for(int j=0;j<n;j++){
			aa[i][j]=a[n*i+j];
		}
	}
   double l[n][n];
   double u[n][n];
   for(int i=1;i<=n;i++){
   	u[0][i-1]=aa[0][i-1];
   	l[i-1][0]=aa[0][i-1]/u[0][0];
   }
   for (int r=2;r<=n;r++){
   	  for(int i=r;i<=n;i++){
   	  	double sum1=0;
   	  	for(int k=1;k<=(r-1);k++){
   	  		sum1=sum1+l[r-1][k-1]*u[k-1][i-1];
			 }
		u[r-1][i-1]=aa[r-1][i-1]-sum1;
   	  	if(r!=n){
   	  		double sum2=0;
   	  		for(int k=1;k<=(r-1);k++){
   	  		sum2=sum2+l[i-1][k-1]*u[k-1][r-1];
			 }
			 l[i-1][r-1]=(aa[i-1][r-1]-sum2)/u[r-1][r-1];
			 }
		 }
   }
   double y[n];
   y[0]=b[0];
   for(int i=2;i<=n;i++){
   	 double sum1=0;
   	 for(int k=1;k<=(i-1);k++){
   	 	sum1=l[i-1][k-1]*y[k-1];
		}
		y[i-1]=b[i-1]-sum1;
   }
   double x[n];
   x[n-1]=y[n-1]/u[n-1][n-1];
   for(int i=(n-1);i>=1;i--){
   	  double sum2=0;
   	  for(int k=(i+1);k<=n;k++){
   	  	 sum2=sum2+u[i-1][k-1]*x[k-1];
		 }
		 x[i-1]=(y[i-1]-sum2)/u[i-1][i-1];
   }
   cout<<"输出求解的U\n";
   for(int i=0;i<n;i++){
   	for(int j=0;j<n;j++){
   		cout<<u[i][j]<<"   ";
	   }
	   cout<<"\n";
   } 
   cout<<"输出求解的L\n";
   for(int i=0;i<n;i++){
   	for(int j=0;j<n;j++){
   		cout<<l[i][j]<<"   ";
	   }
	   cout<<"\n";
   } 
   cout<<"输出求解的x\n";
   for(int i=0;i<n;i++){
   	cout<<x[i]<<"   ";
   }
	return;
}
int main(){
	int n;
	cout<<"输入线性方程组的阶数";
	cin>>n;
	int c=n*n;
	double a[c];
	for(int i=0;i<c;i++){
		cout<<"按行输入矩阵的元素";
		cin>>a[i];
	}
	double b[n];
	for(int i=0;i<n;i++){
		cout<<"输入b的元素";
		cin>>b[i];
	}
	fenjie(a,b,n);
	return 0;
}