lu.cpp

用Doolittle式的LU分解来求解线性方程组。

#include<math.h>
#include<iostream.h>
#include<fstream.h>
//using namespace std;
void main(void)
{
	int LU(double **A, int n);
	int cholesky( double **A, int n );
	int i, j;//循环指标
	
	int n = 0;
	double** A;//矩阵
	
/*	while ( n <= 0 ){
			cout << "Enter the scale of the matrix:" << endl;
			cin >> n;
	}
*/  ///  键盘输入


	
	
//--------------------------//
	ifstream fin ( "input.txt" );
	fin >> n; 
//---------------------------//文件输入
	
	A = new double* [n];
	for ( i = 0; i < n; i++ )
		A[i] = new double[n];//生成矩阵



//------------------------//
	for ( i = 0; i < n; i++ )
		for ( j = 0; j < n; j++ )
			fin >> A[i][j];			//输入矩阵
//----------------------------------//文件输入
                                  

/*	cout << "Enter the Matrix:" << endl;
	for ( i = 0; i < n; i++ )
		for ( j = 0; j < n; j++ )
			cin >> A[i][j];			//输入矩阵
*/   ///  键盘输入
/*
	for ( j = 1; j < n; j++ )
		A[j][0] = A[j][0] / A[0][0];	//计算L的第0列

	for ( i = 1; i < n-1; i++ )
	{
		for ( j = i; j < n; j++ )
			for ( k = 0; k < i; k++ )
				A[i][j] = A[i][j] - A[i][k]*A[k][j]; //计算U的第i列
	    
		for ( j = i+1; j < n; j++ )
		{
			for ( k = 0; k < i; k++ )
				A[j][i] = A[j][i] -A[j][k]*A[k][i];
			if ( fabs(A[i][i]) < 1e-6 )	return 0;
			A[j][i] = A[j][i] / A[i][i];
		}		//计算L的第i列
	}

	for ( k = 0; k < n-1; k++ )
		A[n-1][n-1] = A[n-1][n-1] - A[n-1][k]*A[k][n-1];
*/

//----------------------//
	if( LU(A, n) ) {

	ofstream fout ( "output.txt" );

	fout.precision ( 10 ); // 输出精度
	fout.setf ( ios :: left ); //输出格式
	
		for ( i = 0; i < n; i++ )
		{
			for ( j = 0; j < n; j++ )
			{
				fout.width ( 10 ); //输出格式
				fout << A[i][j];
			}
			fout << endl;
		}
	}
//-----------------------// 文件输出

/*	
	cout.precision ( 10 ); // 输出精度
	cout.setf ( ios :: left ); //输出格式

	for ( i = 0; i < n; i++ )
	{
		for ( j = 0; j < n; j++ )
		{
			cout.width ( 10 ); //输出格式
			cout  << A[i][j];
		}
		cout << endl;
	}
*/ ////屏幕输出
}
int LU(double **A, int n)
{
	int i,j,k;
	for ( j = 1; j < n; j++ )
	A[j][0] = A[j][0] / A[0][0];	//计算L的第0列

	for ( i = 1; i < n-1; i++ )
	{
		for ( j = i; j < n; j++ )
			for ( k = 0; k < i; k++ )
				A[i][j] = A[i][j] - A[i][k]*A[k][j]; //计算U的第i列
	    
		for ( j = i+1; j < n; j++ )
		{
			for ( k = 0; k < i; k++ )
				A[j][i] = A[j][i] -A[j][k]*A[k][i];
			
			if ( fabs(A[i][i]) < 1e-6 )	return 0;
		
			A[j][i] = A[j][i] / A[i][i];
		}		//计算L的第i列
	}

	for ( k = 0; k < n-1; k++ )
		A[n-1][n-1] = A[n-1][n-1] - A[n-1][k]*A[k][n-1];

	return 1;
}