Шпоры по МПиПА / Численные методы / LU-разложение матрицы СЛАУ / Исходник / SLAUGauss

#include "stdio.h"
#include "iostream.h"


// возвращает true если LU-разложение для матрицы,
// состоящей из строк и столбцов с currRowAndColumn п numberOfEquation, было найдено
bool getLUDecomposition( double **matrixU, double **matrixL, int currRowAndColumn, int numberOfEquation ) {
	bool result;
	int i, k;			// переменные циклов 
	double tempItem;	// Вспомогательная переменная

	// если текущий элемент на диагонали равен 0 - LU-разложение не найдено
	result = ( matrixU[currRowAndColumn][currRowAndColumn] != 0 );		
	if ( result && currRowAndColumn < numberOfEquation - 1 ) {
		// определение новых значений элементов матрицы U
		// и получение нового столбца в матрице L
		for ( i = currRowAndColumn + 1; i < numberOfEquation; i++ ) {
			matrixL[i][currRowAndColumn] = matrixU[i][currRowAndColumn] / matrixU[currRowAndColumn][currRowAndColumn];
			tempItem = - matrixL[i][currRowAndColumn];
			for ( k = currRowAndColumn; k < numberOfEquation; k++ ) {
				matrixU[i][k] = matrixU[i][k] + matrixU[currRowAndColumn][k]*tempItem;
			}
		}

		// рекурсивный вызов
		result = getLUDecomposition( matrixU, matrixL, currRowAndColumn + 1, numberOfEquation );
	}

	return result;
}

void main() {
	int i, j;
	int size;
	double **matrixU, **matrixL;

    cout << "Gauss'es method of LU.\nEnter system dimension: ";
    cin >> size;

	matrixU = new double*[size];
	matrixL = new double*[size];
	for ( i = 0; i < size; i++ ) {
		matrixU[i] = new double[size];
		matrixL[i] = new double[size];
		for ( j = 0; j < size; j ++ ){
			matrixL[i][j] = ( i == j ? 1 : 0 );
		}
	}

	for ( i = 0; i < size; i ++ ){
		cout << "Enter " << i + 1 << " row: ";
		for ( j = 0; j < size; j ++ ){
			cin >> matrixU[i][j];
		}
	}

	if ( !getLUDecomposition( matrixU, matrixL, 0, size ) ) {
		cout << "LU-decomposition for this matrix not found";
	}
	else {
		cout << "L-matrix is:\n";
		for ( i = 0; i < size; i ++ ){
			for ( j = 0; j < size; j ++ ){
				cout << matrixL[i][j] << " ";
			}
			cout << "\n";
		}
		cout << "\nU-matrix is:\n";
		for ( i = 0; i < size; i ++ ){
			for ( j = 0; j < size; j ++ ){
				cout << matrixU[i][j] << " ";
			}
			cout << "\n";
		}
	}

    cout << "\nPress \"Enter\" to continue..." << endl; 
    getchar();	
}