void Matrix::eliminate()
{
	/*

	[]-> [][][][]
	[]-> [][][][]
	[]-> [][][][]

	*/

	/*
			Прямой ход
		Задача: Получить верхнюю треугольную матрицу и единицы на главной диагонали
	*/

	std::vector<double> row;

	int totalRows = _size + 1;
	int totalColoumns = _size;

	#pragma omp parallel for
	for (int n = 0; n < _size; n++) {

		/* 1. Search the row with one non-zero value */

		for (int j = n; j < totalRows; j++) {
			for (int i = n; i < totalColoumns; i++) {
				row.push_back(_matrix[i][j]);
			}

			if (std::any_of(row.cbegin(), row.cend(), [](auto a) -> bool { return a != 0; })) {
				if (row[0] == 0) {
					for (int i = n + 1; i < totalColoumns; i++) {
						if (_matrix[i][j] != 0) {
							std::swap(_matrix[i], _matrix[n]);
							std::swap(row[i - n], row[0]);
							break;
						}
					}
				}
				break;
			}
			else {
				row.clear();
				continue;
			}
		}


		/* 2. Divide the each element of first n-coloumn by frist element of choosen row  */
		// Make the 1 on base diagonal

		for (int j = n; j < totalRows; j++) {
			if (!row.empty()) {
				_matrix[n][j] /= row[0]; // TODO: Check for zero value
			}
		}

		/* 3.  */
		// Make the 0 under base diagonal

		for (int i = (n + 1); i < totalColoumns; i++) {
			double firstEl = _matrix[i][n];
			for (int j = n; j < totalRows; j++) {
				_matrix[i][j] -= (_matrix[n][j] * firstEl);
			}
		}

		row.clear();
	}

	/*
			Обратный ход
		Задача: Получить единичную матрицу на месте основной
	*/

	double c = 0; // Коэффицент

	for (int i = _size - 1; i > 0; i--) { // Элемент диагонали
		for (int k = i - 1; k >= 0; k--) { // Строка
			c = _matrix[k][i];
			_matrix[k][i] -= _matrix[i][i] * c;
			_matrix[k][_size] -= _matrix[i][_size] * c;
		}
	}

}