OpenMP – распараллеливание методов численного интегрирования

Главная Форумы Программирование Параллельное программирование Использование OpenMP OpenMP – распараллеливание методов численного интегрирования

В этой теме 0 ответов, 1 участник, последнее обновление  Васильев Владимир Сергеевич 6 мес. назад.

  • Автор
    Сообщения
  • #3930

    Ранее на форуме мы разбирали алгоритм метода прямоугольников и его реализацию на языке С++. Теперь посмотрим как его можно распараллелить с помощью OpenMP (рассмотрим несколько вариантов). Все эти варианты применимы в равной степени ко всем методам численного интегрирования – достаточно заменить функцию на метод трапеций или метод Симпсона.

    Наша программа на С++ содержала два цикла. Один из них обеспечивает достижение заданной точности (количество итераций в нем зависит от особенностей интегрируемой функции, ширины интервала и т.д.), а второй – суммирует N площадей прямоугольников. С помощью конструкции параллельного цикла OpenMP можно распараллелить второй цикл:

    double rectangle_integral(pointFunc f, double a, double b, int n) {
      const double h = (b - a) / n;
      double sum = 0.0;
      #pragma omp parallel for reduction(+:sum)
      for (int i = 0; i < n; i++) {
        sum += f(a + i*h);
      }
      return (sum * h);
    }

    Это наиболее очевидный вариант распараллеливания, однако он имеет недостатки. Потоки создаются в начале параллельной области (#pragma omp parallel) и уничтожаются в ее конце. Создание и уничтожение потоков сопряжено со значительными издержками (ведь это полноценные потоки операционной системы). Эти издержки будут возникать всякий раз, когда основной процесс будет доходить до параллельной области, но ведь количество итераций внешнего цикла заранее не известно – теоретически, эти издержки могут быть существенными.

    Возможен другой подход к распараллеливанию. При этом потоки создаются только один раз – до всех циклов. Затем, между потоками происходит разделение частей одной итерации внешнего цикла (интервал интегрирования делится между потоками поровну). Код внешнего цикла теперь выглядит так:
    bool work = true;

      #pragma omp parallel // 1
      {
        const int thread_count = omp_get_num_threads(),  // 2
            thread_id = omp_get_thread_num();
    
        double thread_h = (b-a)/thread_count;
    
        double loc_a = a+thread_h*thread_id,
            loc_b = loc_a + thread_h;
    
        while (true == work) { // 3
          #pragma omp master  // 4
          {
            s = s1;
            s1 = 0;
            n*=2;
          }
          #pragma omp barrier // 5
          double loc_s = rectangle_integral(f, loc_a, loc_b, ceil(n/thread_count)); // 6
          #pragma omp critical // 7
          {
            s1 += loc_s;
          }
          #pragma omp barrier // 8
          #pragma omp master // 9
          {
            if (fabs(s1 - s) < eps) {
              work = false;
            }
          }
          #pragma omp barrier // 10
        }
      }

    Наиболее интересные фрагменты кода я пронумеровал:

    1. создаются потоки;
    2. один раз каждый поток получает параметры системы – число запущенных потоков, и свой номер. С их помощью вычисляется ширина интервала интегрирования, обрабатываемая одним потоком и границы интервала для текущего потока;
    3. цикл выполняется до тех пор, пока общая для всех потоков переменная work не станет равной false;
    4. главный поток (master) выполняет инициализацию переменной-аккумулятора суммы и рассчитывает новое значение n, отвечающей за количество разбиений интервала интегрирования;
    5. выполняется явная барьерная синхронизация, т.к. директива master не задает неявную синхронизацию, а остальные потоки должны выполнять дальнейшие вычисления только после того, как значение n будет изменено главным потоком;
    6. каждый поток вызывает функцию расчета интеграла на своем интервале (внутри функции никакие потоки теперь не создаются). Для чистоты эксперимента параметр n делится на количество потоков – чтобы количество прямоугольников на всем интервале было равно n, как и в последовательной программе. В результате каждый поток вычисляет свою часть интеграла;
    7. у каждого потока есть своя часть результата, но нам надо сложить их всех в переменную s1, которая является общей. Для этого используем критическую секцию (т.к. она единственная, то пусть она остается безымянной);
    8. поток, вышедший из критической секции продолжает вычисления не дожидаясь пока оттуда выйдет последний поток, иными словами, critical не задает неявной барьерной синхронизации, поэтому нужна явная (с помощь barrier);
    9. кто-то (пусть главный поток, т.е. master) должен решить нужна ли еще одна итерация – если точность достигнута, переменная work получает значение false;
    10. до тех пор, пока главный поток не решит нужно ли продолжать вычисления, другие потоки не должны к ним приступать, поэтому опять используем явную барьерную синхронизацию.

    Время работы разных версий программы сведены в таблицу, при этом запуски проводились для двух функций:

    double f(double x) {
      return x*x+10*x;
    }
    
    double g(double x) {
      return 10*x;
    }

    [0, 100], e=0.005 [0, 200], e=0.005
    последовательная – f 3.54728 27.3215
    параллельный цикл – f 1.85085 13.9019
    ручное распараллеливание – f 1.8132 14.8641
    последовательная – g 0.381472 1.55636
    параллельный цикл – g 0.213399 1.22071
    ручное распараллеливание – g 0.212261 1.30508

    При одних и тех же параметрах вычислительной среды (Процессор: Intel® Pentium(R) CPU G2130 @ 3.20GHz × 2, ОЗУ: 5.8 Гб) время работы программы для функций f и g сильно отличается. Дело не в вычислительной сложности самих функций, а в числе итераций внешнего цикла – чем быстрее растет функция, тем больше разбиений придется выполнить, используя метод прямоугольников. При этом, время работы обоих вариантов распараллеливания отличается мало – если первый вариант имел издержки на создание потоков, то второй сопряжен с примерно такими же по времени издержками на их синхронизацию.

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