矩阵LU的分解

#include<iostream>
#include<iomanip>
#include<fstream>
#include<string>
#include<stdlib.h>
using namespace std;
void main()
{
    //ifstream inf("input.dat");
	int A_N;       //矩阵A的阶数
	cout<<"请输入矩阵的阶数"<<endl;
	cin>>A_N;
	double **A_Matrix;    //矩阵A
	double **L_Matrix;    //矩阵L
	double **U_Matrix;    //矩阵U
	//cout<<"请输入矩阵A的行数和列数"<<endl;
    //cout<<"阶数："<<'	';    
	//A_N=4;         //输入行数
	A_Matrix=(double **)new double *[A_N];// 为矩阵A的行数动态分配空间
	L_Matrix=(double **)new double *[A_N]; //为矩阵L的行数动态分配空间
	U_Matrix=(double **)new double *[A_N]; //为矩阵U的行数动态分配空间
	for(int i=0;i<A_N;i++)
	{
		A_Matrix[i]=new double[A_N];   //为矩阵A的每一列动态分配空间
		L_Matrix[i]=new double[A_N];   //为矩阵L的每一列动态分配空间
		U_Matrix[i]=new double[A_N];   //为矩阵U的每一列动态分配空间
	}
	cout<<"请按行优先的顺序输入矩阵A"<<endl;
	for(int i=0;i<A_N;i++)
	{
		for(int j=0;j<A_N;j++)
		{
			cin>>A_Matrix[i][j];   //矩阵输入
		}
	}
	//cout<<"请以行优先的顺序输入"<<endl;
//	for( i=0;i<A_N;i++)
//	{
        //cout<<"请输入第"<<i<<"行"<<endl;
		//for(int j=0;j<A_N;j++)
	/* int i=0;	//{
	while(!inf.eof())
	{
		char *s1,*s2,*s3,*s4;//声明字符
		s1=new char[4];    //动态分配内存
		s2=new char[4];
		s3=new char[4];
		s4=new char[4];
		s5=new char[2];
        inf>>s1>>s2>>s3>>s4;//将文件input.dat中的数据存入字符
        A_Matrix[i][0]=atof(s1);//将读取到的字符转换为double类型
	    A_Matrix[i][1]=atof(s2);
	    A_Matrix[i][2]=atof(s3);
	    A_Matrix[i][3]=atof(s4);
	    A_Matrix[i][4]=atof(s5);
	    i++;
	}
	inf.close();*/
	   //输入矩阵A
		//}
//	}
	for(int i=0;i<A_N;i++)
	{
		U_Matrix[0][i]=A_Matrix[0][i]; //容易看出A和U的第一行元素对应相等
		for(int j=i+1;j<A_N;j++)     //U矩阵的下三角为0；
		{
			U_Matrix[j][i]=0;	
		}	
	}
	for(int i=0;i<A_N;i++)
	{
		L_Matrix[i][0]=A_Matrix[i][0]/U_Matrix[0][0];//矩阵L为单位下三角矩阵
		                                             //且L的第一列为矩阵A的第一列除以U的第一个元素
		L_Matrix[i][i]=1;         //主对角线为1；
		for(int j=i+1;j<A_N;j++)  //L矩阵的上三角为0
		{
			L_Matrix[i][j]=0;
		}
	}
	for(int k=1;k<A_N;k++)  //推理：已知U的前k-1行和L的前k-1列求U的第k行元素和L的第k列元素
	{
		for(int j=k,i=k+1;j<A_N;j++,i++)
		{
            double U_Sum=0;
			double L_Sum=0;   
			for(int s=0;s<=k-1;s++)
			{
				U_Sum+=L_Matrix[k][s]*U_Matrix[s][j];
			}
			U_Matrix[k][j]=A_Matrix[k][j]-U_Sum;
			if(j==A_N-1)break;
			for(int s=0;s<=k-1;s++)
			{
				L_Sum+=L_Matrix[i][s]*U_Matrix[s][k];
			}
			L_Matrix[i][k]=(A_Matrix[i][k]-L_Sum)/U_Matrix[k][k];
		}	
	}
	for(int i=0;i<A_N;i++)   //输出L矩阵
	{
		for(int j=0;j<A_N;j++)
		{
			cout<<std::left<<setw(10)<<L_Matrix[i][j]<<'	';
		}
         cout<<'
';
	}
	cout<<'
';
	cout<<'
';
	for(int i=0;i<A_N;i++)  //输出U矩阵
	{
		for(int j=0;j<A_N;j++)
		{
			cout<<std::left<<setw(10)<<U_Matrix[i][j]<<'	';
		}
		cout<<'
';
	}
	system("pause");
}