#include <iostream>
#include <vector>
#include <iomanip>
#include <fstream>
#include <cmath>
#include <chrono>

using namespace std;

double computeError(const vector<double>& x, const vector<double>& x_exact);

// Функция для вывода матрицы
void printMatrix(const vector<vector<double>>& matrix) {
    for (const auto& row : matrix) {
        for (double val : row)
            cout << setw(10) << val << " ";
        cout << endl;
    }
}

// Функция для вывода вектора
void printVector(const vector<double>& vec) {
    for (double v : vec)
        cout << setw(10) << v << " ";
    cout << endl;
}

// Чтение матрицы из файла
vector<vector<double>> readMatrixFromFile(const string& filename) {
    ifstream file(filename);
    vector<vector<double>> matrix;
    if (file.is_open()) {
        int n;
        file >> n; // Читаем размерность матрицы
        matrix.resize(n, vector<double>(n));
        for (int i = 0; i < n; i++)
            for (int j = 0; j < n; j++)
                file >> matrix[i][j];
        file.close();
    }
    return matrix;
}

// Чтение вектора из файла
vector<double> readVectorFromFile(const string& filename) {
    ifstream file(filename);
    vector<double> vec;
    if (file.is_open()) {
        int n;
        file >> n; // Читаем размерность вектора
        vec.resize(n);
        for (int i = 0; i < n; i++)
            file >> vec[i];
        file.close();
    }
    return vec;
}

// LU-разложение
bool LU_SQ_Decomposition(const vector<vector<double>>& A, vector<vector<double>>& L, vector<vector<double>>& U) {
    int n = A.size();
    L = vector<vector<double>>(n, vector<double>(n, 0));
    U = vector<vector<double>>(n, vector<double>(n, 0));

    for (int i = 0; i < n; i++) {
        for (int j = 0; j < n; j++) {
            if (j >= i) { // Upper
                U[i][j] = A[i][j];
                for (int k = 0; k < i; k++)
                    U[i][j] -= L[i][k] * U[k][j];
            } else { // Lower
                L[i][j] = A[i][j];
                for (int k = 0; k < j; k++)
                    L[i][j] -= L[i][k] * U[k][j];
                L[i][j] /= U[j][j];
            }
        }
        L[i][i] = 1; // Diagonal elements of L are 1
    }
    return true;
}

// Прямой ход
vector<double> forwardSubstitution(const vector<vector<double>>& L, const vector<double>& F) {
    int n = L.size();
    vector<double> y(n);
    for (int i = 0; i < n; i++) {
        double sum = 0;
        for (int j = 0; j < i; j++)
            sum += L[i][j] * y[j];
        y[i] = (F[i] - sum) / L[i][i];
    }
    return y;
}

// Обратный ход
vector<double> backwardSubstitution(const vector<vector<double>>& U, const vector<double>& y) {
    int n = U.size();
    vector<double> x(n);
    for (int i = n - 1; i >= 0; i--) {
        double sum = 0;
        for (int j = i + 1; j < n; j++)
            sum += U[i][j] * x[j];
        x[i] = (y[i] - sum) / U[i][i];
    }
    return x;
}

// Метод Гаусса
bool gaussianElimination(vector<vector<double>>& A, vector<double>& b) {
    int n = A.size();
    for (int i = 0; i < n; i++) {
        int maxRow = i;
        for (int k = i + 1; k < n; k++) {
            if (abs(A[k][i]) > abs(A[maxRow][i]))
                maxRow = k;
        }
        swap(A[i], A[maxRow]);
        swap(b[i], b[maxRow]);

        for (int k = i + 1; k < n; k++) {
            double factor = A[k][i] / A[i][i];
            for (int j = i; j < n; j++)
                A[k][j] -= factor * A[i][j];
            b[k] -= factor * b[i];
        }
    }

    vector<double> x(n);
    for (int i = n - 1; i >= 0; i--) {
        double sum = 0;
        for (int j = i + 1; j < n; j++)
            sum += A[i][j] * x[j];
        x[i] = (b[i] - sum) / A[i][i];
    }
    // Распечатаем результат здесь для отладки
    cout << "Решение (метод Гаусса):" << endl;
    printVector(x);
    return true;
}

int main() {
    setlocale(LC_ALL, "Russian");

    // Чтение матрицы и вектора из файлов
    vector<vector<double>> A = readMatrixFromFile("matrix.txt");
    vector<double> F = readVectorFromFile("vector.txt");
    int n = A.size();
    vector<vector<double>> L, U;
    vector<double> x_exact(n, 1); // Точное решение

    // Решение с помощью LU-разложения
    cout << "Решение для матрицы из файла matrix.txt:" << endl;
    auto start = chrono::high_resolution_clock::now();
    if (LU_SQ_Decomposition(A, L, U)) {
        vector<double> y = forwardSubstitution(L, F);
        vector<double> x = backwardSubstitution(U, y);
        double error = computeError(x, x_exact);
        cout << "Матрица A:" << endl;
        printMatrix(A);
        cout << "Решение LU-разложения:" << endl;
        printVector(x);
        cout << "Погрешность решения: " << error << endl;
    }
    auto end = chrono::high_resolution_clock::now();
    cout << "Время выполнения LU-разложения: " << chrono::duration<double>(end - start).count() << " секунд" << endl;

    // Метод Гаусса
    vector<vector<double>> A_gauss = A;
    vector<double> F_gauss = F;
    cout << "\nРешение с использованием метода Гаусса:" << endl;
    start = chrono::high_resolution_clock::now();
    if (gaussianElimination(A_gauss, F_gauss)) {
        cout << "Метод Гаусса успешно выполнен. Решение:" << endl;
        printVector(F_gauss);
    }
    end = chrono::high_resolution_clock::now();
    cout << "Время выполнения метода Гаусса: " << chrono::duration<double>(end - start).count() << " секунд" << endl;

    return 0;
}

fix this code