Шпоры по МПиПА / Численные методы / Поиск обратной матрицы методом Гаусса / Исходник / SLAUGauss

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

// приведение матрицы коэффициентов к диагональноиу виду и формирование b-матрицы
// для дальнейшего использования ее для нахождения обратной матрицы
bool prepareBMatrix( double **coefficients, double **bMatrix, int currRowAndColumn, int numberOfEquation ) {
	bool result;
	bool allElementsInCurrentColumnEqualsZero = true;
	int i, k, row;
	double tempItem;

	if ( currRowAndColumn == numberOfEquation - 1 ) {
		result = ( coefficients[currRowAndColumn][currRowAndColumn] != 0 );
	}
	else {
		// если в текущем столбце все элементы ниже текущей строки равны 0  - решение не может быть найдено
		// если элемент на диагонали - нулевой, ищем в текущем столбце ненулевой элемент в строках ниже текущей
		// для перестановки
		for ( i = currRowAndColumn; i < numberOfEquation; i++ ) {
			if ( coefficients[i][currRowAndColumn] != 0 ) {
				allElementsInCurrentColumnEqualsZero = false;
				row = i;
				break;
			}
		}

		if ( allElementsInCurrentColumnEqualsZero ) {
			result = false;
		}
		else {

			// если нужно - переставляем местами строки в матрицах 
			if ( row != currRowAndColumn ) {
				for ( i = currRowAndColumn; i < numberOfEquation; i++ ) {
					tempItem = coefficients[currRowAndColumn][i];
					coefficients[currRowAndColumn][i] = coefficients[row][i];
					coefficients[row][i] = tempItem;
				}
				for ( i = 0; i < numberOfEquation; i++ ) {
					tempItem = bMatrix[currRowAndColumn][i];
					bMatrix[currRowAndColumn][i] = bMatrix[row][i];
					bMatrix[row][i] = tempItem;
				}
			}

			// пересчитываем элементы в матрицах
			for ( i = currRowAndColumn + 1; i < numberOfEquation; i++ ) {
				tempItem = -coefficients[i][currRowAndColumn] / coefficients[currRowAndColumn][currRowAndColumn];
				for ( k = currRowAndColumn; k < numberOfEquation; k++ ) {
					coefficients[i][k] = coefficients[i][k] + coefficients[currRowAndColumn][k]*tempItem;
				}
				for ( k = 0; k < numberOfEquation; k++ ) {
					bMatrix[i][k] = bMatrix[i][k] + bMatrix[currRowAndColumn][k]*tempItem;
				}
			}

			// рекурсивный вызов с увеличенным счетчиком текущих столбца и строки
			result = prepareBMatrix( coefficients, bMatrix, currRowAndColumn + 1, numberOfEquation );
		}
	}

	return result;
}

// было ли найдено решение, если да - итог в параметре solution
bool inverseMatrix( double **coefficients, int numberOfEquation, double **solution ) {
	bool result;
	int i, k, j;
	double **bMatrix = new double*[numberOfEquation];
	for ( i = 0; i < numberOfEquation; i ++ ){
		bMatrix[i] = new double[numberOfEquation];
		for ( k = 0; k < numberOfEquation; k ++ ){
			bMatrix[i][k] = ( i == k ? 1 : 0 );
		}
	}
	
	// преобразование входной матрицы в диагональную и формирование вспомогательной b-матрицы
	result = prepareBMatrix( coefficients, bMatrix, 0, numberOfEquation );

	// собственно нахождение обратной матрицы из numberOfEquation систем по numberOfEquation уравнений в каждой
	if ( result ) {
		for ( j = 0; j < numberOfEquation; j ++ ) {
			for ( i = numberOfEquation - 1; i >= 0; i -- ) {
				solution[i][j] = bMatrix[i][j];
				for ( k = i + 1; k < numberOfEquation; k ++ ) {
					solution[i][j] = solution[i][j] - solution[k][j]*coefficients[i][k];
				}
				solution[i][j] = solution[i][j] / coefficients[i][i];
			}
		}
	}

	return result;
}


void main() {
	int i, j;
	int size;
	double **coefficients, **solution;

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

	coefficients = new double*[size];
	solution = new double*[size];
	for ( i = 0; i < size; i++ ) {
		coefficients[i] = new double[size];
		solution[i] = new double[size];
	}

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

	if ( !inverseMatrix( coefficients, size, solution ) ) {
		cout << "Solution for this matrix of coefficients not exist";
	}
	else {
		cout << "Inverse matrix is:\n";
		for ( j = 0; j < size; j ++ ){
			for ( i = 0; i < size; i ++ ){
				cout << solution[j][i] << " ";
			}
			cout << "\n";
		}
	}

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