#include <fstream>
#include <iostream>
#include <omp.h>
#include <windows.h>
#include <complex>
using namespace std;

void jacobi(int N) {
    const double eps = 0.0001;
    char nameA[50], nameB[50], matrx[50];
    sprintf_s(nameA, "matrixAComplex%d.txt", N);
    sprintf_s(nameB, "matrixBComplex%d.txt", N);
    sprintf_s(matrx, "JacobiRes%d.txt", N);
    ifstream matrA(nameA);
    ifstream matrB(nameB);
    ofstream matrX(matrx);
    if (!matrA.is_open() || !matrB.is_open() || !matrX.is_open()) {
        cout << "Ошибка открытия файла" << endl;
        return;
    }
    complex<double>* matrixB = new complex<double>[N];
    complex<double>* matrixx = new complex<double>[N];
    complex<double>* matrixxx = new complex<double>[N];
    complex<double>** matrixA = new complex<double>* [N];

    for (int i = 0; i < N; i++) {
        matrixA[i] = new complex<double>[N];
    }

    for (int i = 0; i < N; i++) {
        for (int j = 0; j < N; j++) {
            matrA >> matrixA[i][j];
        }
        matrB >> matrixB[i];
        matrixx[i] = 1.0; // Инициализация начального приближения
        matrixxx[i] = 1.0; // Инициализация вектора для нового приближения
    }

    matrA.close();
    matrB.close();
    double x_max = -1.0, count;
    complex<double> sum1, sum2;
    double tn = omp_get_wtime();
    do {
        x_max = -1.0;
        for (int i = 0; i < N; i++) {
            sum1 = 0.0; // Сумма для элементов ниже диагонали
            sum2 = 0.0; // Сумма для элементов выше диагонали
            for (int j = 0; j < N; j++)
            {
                if (j < i)
                    sum1 += matrixA[i][j] * matrixx[j];
                if (j > i)
                    sum2 += matrixA[i][j] * matrixx[j];
            }
            matrixxx[i] = 1.0 / matrixA[i][i] * (matrixB[i] - sum1 - sum2); // Вычисление нового приближения
            count = abs(matrixxx[i] - matrixx[i]); // Вычисление максимального изменения компонента вектора
            if (x_max < count) {
                x_max = count;
            }
        }
        // Обновление вектора приближений
        for (int i = 0; i < N; i++) {
            matrixx[i] = matrixxx[i];
        }
    } while (x_max > eps);
    double tk = omp_get_wtime();
    cout << "Время выполнения " << N << ": " << tk - tn << endl;
    for (int i = 0; i < N; i++) {
        matrX << matrixx[i] << endl;
    }
    matrA.close();
    matrB.close();
    matrX.close();
    for (int i = 0; i < N; i++)
    {
        delete[] matrixA[i];

    }
    delete[] matrixA;
    delete[] matrixB;
    delete[] matrixx;
    delete[] matrixxx;
}

void gausZeyd(int n) {
    const double eps = 0.0001;
    char nameA[50], nameB[50], matrx[50];
    sprintf_s(nameA, "matrixAComplex%d.txt", n);
    sprintf_s(nameB, "matrixBComplex%d.txt", n);
    sprintf_s(matrx, "GausZeydRes%d.txt", n);
    ifstream matrA(nameA);
    ifstream matrB(nameB);
    ofstream matrX(matrx);
    if (!matrA.is_open() || !matrB.is_open() || !matrX.is_open()) {
        cout << "Ошибка открытия файла" << endl;
        return;
    }

    complex<double>** matrixA = new complex<double>* [n];
    complex<double>* matrixB = new complex<double>[n];
    complex<double>* matrixX = new complex<double>[n];

    for (int i = 0; i < n; i++) matrixA[i] = new complex<double>[n];

    for (int i = 0; i < n; i++)
    {
        for (int j = 0; j < n; j++)
        {
            matrA >> matrixA[i][j];
        }
        matrB >> matrixB[i];
        matrixX[i] = 1.0; // Инициализация начального вектора приближения
    }

    matrA.close();
    matrB.close();

    double xmax;
    complex<double> sum1, sum2;
    double tn = omp_get_wtime();
    double tmp = 0.0;
    complex<double> xv;

    tn = omp_get_wtime();
    do {
        xmax = -1.0;
        for (int i = 0; i < n; i++)
        {
            xv = matrixX[i];
            sum1 = 0.0;
            sum2 = 0.0;
            for (int j = 0; j < n; j++)
            {
                if (j < i)
                    sum1 += matrixA[i][j] * matrixX[j];
                if (j > i)
                    sum2 += matrixA[i][j] * matrixX[j];
            }
            matrixX[i] = (1.0 / matrixA[i][i]) * (matrixB[i] - sum1 - sum2);
            tmp = abs(matrixX[i] - xv);
            if (tmp > xmax) xmax = tmp;
        }
    } while (xmax > eps);

    double tk = omp_get_wtime();
    cout << "Время выполнения " << n << ": " << tk - tn << endl;

    for (int i = 0; i < n; i++) {
        matrX << matrixX[i] << endl;
    }

    matrX.close();

    for (int i = 0; i < n; i++)
    {
        delete[] matrixA[i];

    }
    delete[] matrixA;
    delete[] matrixB;
    delete[] matrixX;
}

void gaus(int n)
{
    const double eps = 0.0001;
    char nameA[50], nameB[50], matrx[50];
    sprintf_s(nameA, "matrixAComplex%d.txt", n);
    sprintf_s(nameB, "matrixBComplex%d.txt", n);
    sprintf_s(matrx, "GausRes%d.txt", n);
    ifstream matrA(nameA);
    ifstream matrB(nameB);
    ofstream matrX(matrx);

    if (!matrA.is_open() || !matrB.is_open() || !matrX.is_open()) {
        cout << "Ошибка открытия файла" << endl;
        return;
    }

    complex<double>** matrixA = new complex<double>* [n];
    complex<double>* matrixB = new complex<double>[n];
    complex<double>* matrixX = new complex<double>[n];

    for (int i = 0; i < n; i++) {
        matrixA[i] = new complex<double>[n];
    }

    for (int i = 0; i < n; i++)
    {
        for (int j = 0; j < n; j++)
        {
            matrA >> matrixA[i][j];
        }
        matrB >> matrixB[i];
        matrixX[i] = 1.0; // Инициализация начального вектора приближения
    }

    matrA.close();
    matrB.close();

    complex<double> w = 0.0, sum = 0.0;
    double tn = omp_get_wtime();

    for (int i = 0; i <= n - 1; i++)
    {
        for (int j = i + 1; j <= n - 1; j++)
        {
            w = matrixA[j][i] / matrixA[i][i];
            for (int k = i; k <= n - 1; k++)
            {
                matrixA[j][k] = matrixA[j][k] - w * matrixA[i][k];
            }
            matrixB[j] = matrixB[j] - w * matrixB[i];
        }
    }

    for (int i = n - 1; i >= 0; i--)
    {
        sum = 0.0;
        for (int j = i; j <= n - 1; j++)
        {
            if (i != NULL)
                sum = sum + matrixA[i][j + 1] * matrixX[j + 1];
        }
        matrixX[i] = (matrixB[i] - sum) / matrixA[i][i];
    }

    double tk = omp_get_wtime();
    cout << "Время выполнения " << n << ": " << tk - tn << endl;

    for (int i = 0; i < n; i++) {
        matrX << matrixX[i] << endl;
    }

    matrX.close();

    for (int i = 0; i < n; i++)
    {
        delete[] matrixA[i];

    }
    delete[] matrixA;
    delete[] matrixB;
    delete[] matrixX;
}

int main()
{
    SetConsoleOutputCP(1251);
    jacobi(100);
    jacobi(250);
    jacobi(500);

    gausZeyd(100);
    gausZeyd(250);
    gausZeyd(500);

    gaus(100);
    gaus(250);
    gaus(500);
}