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

// Выбор точности: раскомментируйте нужное определение
typedef double precision_t;
// typedef float precision_t;

using namespace std;

// Функция для вычисления ошибки между найденным и точным решениями
template <typename T>
double computeError(const vector<T>& x, const vector<T>& x_exact) {
    if (x.size() != x_exact.size()) {
        cerr << "Error: Векторы x и x_exact должны иметь одинаковый размер." << endl;
        return -1.0;
    }
    double error = 0.0;
    for (size_t i = 0; i < x.size(); ++i) {
        error += pow(static_cast<double>(x[i] - x_exact[i]), 2);
    }
    return sqrt(error);
}

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

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

// Чтение матрицы из файла в профильном формате
template <typename T>
vector<vector<T>> readMatrixFromFile(const string& filename) {
    ifstream file(filename);
    vector<vector<T>> matrix;
    if (file.is_open()) {
        int n;
        file >> n; // Читаем размерность матрицы
        matrix.resize(n, vector<T>(n, 0));
        // Чтение матрицы в профильном формате
        for (int i = 0; i < n; ++i) {
            for (int j = 0; j < n; ++j) {
                file >> matrix[i][j];
            }
        }
        file.close();
    } else {
        cerr << "Не удалось открыть файл матрицы: " << filename << endl;
    }
    return matrix;
}

// Чтение вектора из файла
template <typename T>
vector<T> readVectorFromFile(const string& filename) {
    ifstream file(filename);
    vector<T> vec;
    if (file.is_open()) {
        int n;
        file >> n; // Читаем размерность вектора
        vec.resize(n);
        for (int i = 0; i < n; ++i)
            file >> vec[i];
        file.close();
    } else {
        cerr << "Не удалось открыть файл вектора: " << filename << endl;
    }
    return vec;
}

// Запись вектора в файл
template <typename T>
void writeVectorToFile(const string& filename, const vector<T>& vec) {
    ofstream file(filename);
    if (file.is_open()) {
        file << vec.size() << endl;
        for (auto val : vec)
            file << val << " ";
        file << endl;
        file.close();
    } else {
        cerr << "Не удалось открыть файл для записи: " << filename << endl;
    }
}

// LU-разложение с хранением L и U в исходной матрице
template <typename T>
bool LU_Decomposition(vector<vector<T>>& A, vector<int>& pivot, unsigned long long& operation_count) {
    int n = A.size();
    pivot.resize(n);
    for (int i = 0; i < n; ++i)
        pivot[i] = i;

    for (int k = 0; k < n; ++k) {
        // Поиск главного элемента
        T max = abs(A[k][k]);
        int maxRow = k;
        for (int i = k + 1; i < n; ++i) {
            if (abs(A[i][k]) > max) {
                max = abs(A[i][k]);
                maxRow = i;
            }
            operation_count++;
        }

        if (max == 0) {
            cerr << "Матрица вырождена!" << endl;
            return false;
        }

        // Перестановка строк
        if (maxRow != k) {
            swap(A[k], A[maxRow]);
            swap(pivot[k], pivot[maxRow]);
            operation_count += 3 * n; // Предполагаем, что swap требует 3n операций
        }

        // Приведение к верхнетреугольному виду
        for (int i = k + 1; i < n; ++i) {
            A[i][k] /= A[k][k];
            operation_count++;
            for (int j = k + 1; j < n; ++j) {
                A[i][j] -= A[i][k] * A[k][j];
                operation_count += 2;
            }
        }
    }
    return true;
}

// Прямой ход для решения L*y = b с учетом перестановок
template <typename T>
vector<T> forwardSubstitution(const vector<vector<T>>& A, const vector<T>& b, const vector<int>& pivot, unsigned long long& operation_count) {
    int n = A.size();
    vector<T> y(n, 0);

    for (int i = 0; i < n; ++i) {
        y[i] = b[pivot[i]];
        for (int j = 0; j < i; ++j) {
            y[i] -= A[i][j] * y[j];
            operation_count += 2;
        }
        // L имеет единицы на диагонали
    }
    return y;
}

// Обратный ход для решения U*x = y
template <typename T>
vector<T> backwardSubstitution(const vector<vector<T>>& A, const vector<T>& y, unsigned long long& operation_count) {
    int n = A.size();
    vector<T> x(n, 0);

    for (int i = n - 1; i >= 0; --i) {
        x[i] = y[i];
        for (int j = i + 1; j < n; ++j) {
            x[i] -= A[i][j] * x[j];
            operation_count += 2;
        }
        x[i] /= A[i][i];
        operation_count++;
    }
    return x;
}

// Метод Гаусса с выбором главного элемента и подсчетом операций
template <typename T>
bool GaussianElimination(vector<vector<T>> A, vector<T> b, vector<T>& x, unsigned long long& operation_count) {
    int n = A.size();
    vector<int> pivot(n);
    for (int i = 0; i < n; ++i)
        pivot[i] = i;

    // Прямой ход
    for (int k = 0; k < n; ++k) {
        // Выбор главного элемента
        T max = abs(A[k][k]);
        int maxRow = k;
        for (int i = k + 1; i < n; ++i) {
            if (abs(A[i][k]) > max) {
                max = abs(A[i][k]);
                maxRow = i;
            }
            operation_count++;
        }

        if (max == 0) {
            cerr << "Матрица вырождена!" << endl;
            return false;
        }

        // Перестановка строк
        if (maxRow != k) {
            swap(A[k], A[maxRow]);
            swap(b[k], b[maxRow]);
            swap(pivot[k], pivot[maxRow]);
            operation_count += 3 * n; // Предполагаем, что swap требует 3n операций
        }

        // Нули под диагональю
        for (int i = k + 1; i < n; ++i) {
            T factor = A[i][k] / A[k][k];
            operation_count++;
            A[i][k] = 0;
            for (int j = k + 1; j < n; ++j) {
                A[i][j] -= factor * A[k][j];
                operation_count += 2;
            }
            b[i] -= factor * b[k];
            operation_count += 1;
        }
    }

    // Обратный ход
    x.assign(n, 0);
    for (int i = n - 1; i >= 0; --i) {
        if (A[i][i] == 0) {
            cerr << "Нулевой элемент на диагонали в строке " << i << endl;
            return false;
        }
        x[i] = b[i];
        for (int j = i + 1; j < n; ++j) {
            x[i] -= A[i][j] * x[j];
            operation_count += 2;
        }
        x[i] /= A[i][i];
        operation_count++;
    }
    return true;
}

int main() {
    // Файлы ввода и вывода
    string matrixFile = "matrix.txt";
    string vectorFile = "vector.txt";
    string outputFile_LU = "solution_LU.txt";
    string outputFile_Gauss = "solution_Gauss.txt";

    // Чтение матрицы и вектора
    vector<vector<precision_t>> A = readMatrixFromFile<precision_t>(matrixFile);
    vector<precision_t> b = readVectorFromFile<precision_t>(vectorFile);

    if (A.empty() || b.empty()) {
        cerr << "Ошибка при чтении файлов. Проверьте содержимое." << endl;
        return 1;
    }

    int n = A.size();
    if (b.size() != n) {
        cerr << "Размер вектора правой части не соответствует размерности матрицы." << endl;
        return 1;
    }

    // Задание точного решения (предполагается, что оно известно)
    // В этом примере x_exact = (1, 2, ..., n)
    vector<precision_t> x_exact(n, 0);
    for (int i = 0; i < n; ++i)
        x_exact[i] = i + 1;

    // Клонирование матрицы A для LU-разложения
    vector<vector<precision_t>> A_LU = A;
    vector<int> pivot;
    unsigned long long operations_LU = 0;

    // LU-разложение
    cout << "Выполнение LU-разложения..." << endl;
    auto start_LU = chrono::high_resolution_clock::now();
    bool success_LU = LU_Decomposition<precision_t>(A_LU, pivot, operations_LU);
    auto end_LU = chrono::high_resolution_clock::now();
    double time_LU = chrono::duration<double>(end_LU - start_LU).count();

    if (!success_LU) {
        cerr << "LU-разложение не удалось." << endl;
        return 1;
    }

    // Прямой и обратный ход
    unsigned long long operations_forward = 0;
    vector<precision_t> y = forwardSubstitution<precision_t>(A_LU, b, pivot, operations_forward);

    unsigned long long operations_backward = 0;
    vector<precision_t> x_LU = backwardSubstitution<precision_t>(A_LU, y, operations_backward);

    // Вычисление полной погрешности
    double error_LU = computeError<precision_t>(x_LU, x_exact);

    // Запись решения LU-разложения в файл
    writeVectorToFile<precision_t>(outputFile_LU, x_LU);

    // Вывод результатов LU-разложения
    cout << "Результаты LU-разложения:" << endl;
    cout << "Время выполнения: " << fixed << setprecision(6) << time_LU << " секунд" << endl;
    cout << "Количество операций: " << operations_LU + operations_forward + operations_backward << endl;
    cout << "Погрешность: " << fixed << setprecision(6) << error_LU << endl;

    // Метод Гаусса
    vector<vector<precision_t>> A_Gauss = A; // Копия матрицы для метода Гаусса
    vector<precision_t> b_Gauss = b;         // Копия вектора правой части
    vector<precision_t> x_Gauss;
    unsigned long long operations_Gauss = 0;

    cout << "\nВыполнение метода Гаусса..." << endl;
    auto start_Gauss = chrono::high_resolution_clock::now();
    bool success_Gauss = GaussianElimination<precision_t>(A_Gauss, b_Gauss, x_Gauss, operations_Gauss);
    auto end_Gauss = chrono::high_resolution_clock::now();
    double time_Gauss = chrono::duration<double>(end_Gauss - start_Gauss).count();

    if (!success_Gauss) {
        cerr << "Метод Гаусса не удалось." << endl;
        return 1;
    }

    // Вычисление полной погрешности для метода Гаусса
    double error_Gauss = computeError<precision_t>(x_Gauss, x_exact);

    // Запись решения метода Гаусса в файл
    writeVectorToFile<precision_t>(outputFile_Gauss, x_Gauss);

    // Вывод результатов метода Гаусса
    cout << "Результаты метода Гаусса:" << endl;
    cout << "Время выполнения: " << fixed << setprecision(6) << time_Gauss << " секунд" << endl;
    cout << "Количество операций: " << operations_Gauss << endl;
    cout << "Погрешность: " << fixed << setprecision(6) << error_Gauss << endl;

    // Сравнение результатов
    cout << "\nСравнение LU-разложения и метода Гаусса:" << endl;
    cout << "Погрешность LU: " << fixed << setprecision(6) << error_LU << endl;
    cout << "Погрешность Гаусса: " << fixed << setprecision(6) << error_Gauss << endl;
    cout << "Количество операций LU: " << operations_LU + operations_forward + operations_backward << endl;
    cout << "Количество операций Гаусса: " << operations_Gauss << endl;

    return 0;
}

Вот мой код. В каком формате представлена матрица (ленточном или профильном)?
Можешь проанализировать. Выполняется ли тут LU(sq) разложение
и обнови метод гаусса так, чтобы он был с выбором ведущего элемента, для плотных матриц.