Распараллеливание задачи численного интегрирования

Помечено: 

В этой теме 3 ответа, 2 участника, последнее обновление 1 месяц назад.

  • Автор
    Сообщения
  • #5363
    @sss321

    Найти численное решение определенного интеграла методом Симпсона. Промежуток интегрирования разбить на несколько частей, применить на каждой части формулу Симпсона (соответственно каждая часть — отдельный процесс). Провести вычисления не разбивая отрезок, затем разбить на 2 части, потом на 4 и затем на 16. Полученные решения сравнить с точным решением.

  • #5364
    @admin

    Итак, берем готовое решение — Метод парабол (Симпсона) на С++ и распараллеливаем.

    Собственно сделать это можно так:

    double simpson_integral(pointFunc f, double a, double b, int n) {
      const double h = (b-a)/n;
      double k1 = 0, k2 = 0;
      #pragma omp parallel for reduction(+:k1,k2)
      for(int i = 1; i < n; i += 2) {
        k1 += f(a + i*h);
        k2 += f(a + (i+1)*h);
      }
      return h/3*(f(a) + 4*k1 + 2*k2);
    }

    Тут добавилась единственная строка. Как это работает?
    Директива Parallel порождает M потоков, где M по умолчанию равно количеству ядер в вашей системе.
    Итерации цикла (грубо говоря, это и есть отрезки) распределяются между потоками — по умолчанию примерно поровну. Каждый из них считает свое локальное значение k1 и k2, затем все результаты складываются в глобальные переменные с таким же именем. Это должно хорошо работать.

    Чтобы разбить на 1, 2, 4, 16 частей — можно использовать дополнительную опцию num_threads(M):
    #pragma omp parallel for reduction(+:k1,k2) num_threads(4)

    Собственно можно провести такие эксперименты.

  • #5365
    @sss321

    А если сначала сразу разбить отрезок на 16 частей, это значит задавать M=16, чтобы в формулу h=(b-a)/M передалось значение 16?
    И тогда потом не нужно увеличивать число шагов в 2 раза в цикле do…while ?

    • #5366
      @admin

      Чтобы «сразу разбить»:

      double h = (b-a)/M;
      double sum = 0;
      #pragma omp parallel num_threads(M) reduction(+:sum)
      {
        int thread_num = omp_get_thread_num();
        double local_a = a + h*thread_num;
        double local_b = local_a + h;
        sum = simpson_integral(foo, local_a, local_b, n);
      }

      Должно работать. Функция интегрирования при этом используется последовательная. Мы создаем M потоков, делим весь интервал [a, b] на M частей и в каждом потоке обрабатываем свою часть. С помощью опции reduction все полученные результаты складываются в общую переменную.
      Тут нужно учитывать, что интервал теперь будет биться не на n участков, а на n*M.

Для ответа в этой теме необходимо авторизоваться.