Yacobi_OpenMP

#include <iostream>
#include <omp.h>
#include <cmath>
#include <complex>
#define pi 3.14159265
#define eps 0.000001
using namespace std;

complex<double> gamma(double n, const double& K)
{
    complex<double> res;
    if (K * K >= pi * pi * n * n)
    {
        res = sqrt(K * K - pi * pi * n * n);
    }
    else
    {
        res = 1i * sqrt(pi * pi * n * n - K * K);
    }
    return res;
}
complex<double> Jnm(double n, double m, const double& a0, const double& a1)
{
    complex<double> res;
    if (n == m)
    {
        res = a1 - a0 - ((sin((n + m) * pi * a1) - sin((n + m) * pi * a0)) / ((n + m) * pi));
    }
    else
    {
        res = (sin((n - m) * pi * a1) - sin((n - m) * pi * a0)) / ((n - m) * pi) - ((sin((n + m) * pi * a1) - sin((n + m) * pi * a0)) / ((n + m) * pi));
    }
    return res;
}
complex<double> Inm(double n, double m, const double& a0, const double& a1)
{
    complex<double> res;
    if (n == m)
    {
        res = 1.0 - Jnm(n, m, a0, a1);
    }
    else
    {
        res = -Jnm(n, m, a0, a1);
    }
    return res;
}
complex<double> sumM(int i, int j, const int& N, const double& c, const double& a0, const double& a1, const double& K)
{
    complex<double>s;
    for (int l = 0; l < N; l++)
        s += gamma(l + 1, K) * Jnm(j + 1, l + 1, a0, a1) * Inm(l + 1, i + 1, a0, a1) / (1.0 - exp(2.0 * 1i * gamma(l + 1, K) * c));
    return s;
}
void getAandb(const int&n, complex<double>** matrA, complex<double>* vecb, const int& N, const double& l, const double& K, const double& a0, const double& a1, const double& c)
{
    int i, j;
    #pragma omp parallel num_threads(n) private(i,j)
    {
        #pragma omp for schedule(dynamic)
        for (i = 0; i < N; i++)
        {
            for (j = 0; j < N; j++)
            {
                if (i == j)
                {
                    matrA[i][j] = gamma(i + 1, K) - (1.0 - exp(2.0 * 1i * gamma(i + 1, K) * c)) * sumM(i, j, N, c, a0, a1, K);
                }
                else
                {
                    matrA[i][j] = -(1.0 - exp(2.0 * 1i * gamma(j + 1, K) * c)) * sumM(i, j, N, c, a0, a1, K);
                }
            }
        }
        #pragma omp for
        for (int i = 0; i < N; i++)
            vecb[i] = gamma(l, K) * Inm(l, i + 1, a0, a1);
    }
}
void show(complex<double>* x, int N)
{
    for (int i = 0; i < N; i++)
        cout << abs(x[i]) << "\n";
    cout << "\n";
}
void Yacobi(const int& n, const int& N, const double& c, const double& l, const double& a0, const double& a1, const double& K)
{
    complex<double>** matrA = new complex<double>* [N];//матрица коэффициентов
    for (int i = 0; i < N; i++)
        matrA[i] = new complex<double>[N];
    complex<double>* F = new complex<double>[N];//столбец свободных членов
    getAandb(n, matrA, F, N, l, K, a0, a1, c);
    complex<double>* x = new complex<double>[N]; //начальное приближение, ответ также записывается в х
    complex<double>* tempx = new complex<double>[N];
    double norm;
    do
    {
        for (int i = 0; i < N; i++)
        {
            tempx[i] = F[i];
            for (int g = 0; g < N; g++)
            {
                if (i != g)
                    tempx[i] -= matrA[i][g] * x[g];
            }
            tempx[i] /= matrA[i][i];
        }
        norm = abs(x[0] - tempx[0]);
        for (int h = 0; h < N; h++)
        {
            if (abs(x[h] - tempx[h]) > norm)
                norm = abs(x[h] - tempx[h]);
            x[h] = tempx[h];
        }
    } while (norm > eps);
    show(x, N);
    delete[]x;
    delete[]tempx;
    delete[]F;
    for (int i = 0; i < N; i++)
        delete[]matrA[i];
    delete[]matrA;
}
int main()
{
    setlocale(LC_ALL, "RU");
    double c, l, a0, a1, K;
    cout << "Введите параметры с, l, a0, a1, к cоответственно: \n";
    cin >> c >> l >> a0 >> a1 >> K;
    int n, N;
    cout << "Введите размерность системы N: \n";
    cin >> N;
    cout << "Задайте количество потоков: \n";
    cin >> n;
    double start_time, end_time;
    start_time = omp_get_wtime();
    Yacobi(n, N, c, l, a0, a1, K);
    end_time = omp_get_wtime();
    cout << end_time - start_time << "\n";
    return 0;
}