矩阵的lu分解法.cpp

来自「数值分析算法矩阵的LU分解法」· C++ 代码 · 共 127 行

/**************** LU分解法 *******************/
#define AL 100     //局阵最大行
#define ROW 100		//局阵最大列
#include<iostream.h>
#include<stdlib.h>
#include<math.h>
#include<fstream.h>

void main(){
	int n;//阶
	double A[AL][ROW];//系数矩阵
	double b[AL];//右端项
	int i,j;
	char flag='n';//选择标志

	//.......  文件读入  .........
	/*cout<<"是否从文件读入(Y/N):  ";
	cin>>flag;
	if(flag=='y'||flag=='Y'){
		cout<<"请输入文件名:(file.txt)  ";
		char file[10];
		//已存在文件:file1.txt,file2.txt
		cin>>file;  
		cout<<endl;

		ifstream fin(file);

	cout<<"输入阶:  "<<endl;
	fin>>n;
	cout<<n<<endl;

	cout<<"输入系数矩阵A:  "<<endl;
	for(i=0;i<n;i++){
		for(j=0;j<n;j++){
			fin>>A[i][j];
			cout<<A[i][j]<<"  ";
		}
		cout<<endl;
	}
	
	cout<<"输入右端项b:  "<<endl;
	for(i=0;i<n;i++){
		fin>>b[i];
		cout<<b[i]<<"  ";
	}
		cout<<endl;
		cout<<endl;

	}*/


	//........  手工输入  ............
//	else{
	cout<<"输入阶:  "<<endl;
	cin>>n;

	cout<<"输入系数矩阵A:  "<<endl;
	for(i=0;i<n;i++)
		for(j=0;j<n;j++)
			cin>>A[i][j];
	
	cout<<"输入右端项b:  "<<endl;
	for(i=0;i<n;i++)
		cin>>b[i];
	cout<<endl;
//	}
	

//.........计算L矩阵和U矩阵(都放到A矩阵里)..........
	int k;
//	int m;
	double sum_1=0;
	for(k=0;k<n-1;k++){//循环

//L矩阵算法
		for(i=k+1;i<n;i++){
			sum_1=0;
			for(int m=0;m<k;m++)
				sum_1+=A[i][m]*A[m][k];
			A[i][k]=(A[i][k]-sum_1)/A[k][k];
		}
//U矩阵算法
		for(i=k+1;i<n;i++){
			sum_1=0;
			for(int m=0;m<k+1;m++)
				sum_1+=A[k+1][m]*A[m][i];
			A[k+1][i]-=sum_1;
		}

	}
			
//..............  LUX=b  ......................

	//第一步:UX=X_*,算LX_*=b,得到Y(即b[i])
	for(i=0;i<n;i++){
		sum_1=0;
		for(j=i;j>0;j--)
			sum_1+=A[i][j-1]*b[j-1];
		b[i]=b[i]-sum_1;
	}
	//第二步:UY=b,得出答案b[i]
	for(i=n-1;i>=0;i--){
					sum_1=0;
					for(j=i+1;j<n;j++){					
						sum_1=sum_1+A[i][j]*b[j];
				}
					b[i]=(b[i]-sum_1)/A[i][i];
			}

	cout<<"\n答案输出:  ";

	for(i=0;i<n;i++)
		cout<<b[i]<<"  ";
	cout<<endl;

//L矩阵和U矩阵都仍旧存放在A矩阵里
/*	for(i=0;i<n;i++){
		for(j=0;j<n;j++)
			cout<<A[i][j]<<"  ";
		cout<<endl;*/
	
		
}