Распараллеливание алгоритма триангуляции матрицы с OpenMP

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

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

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

    В алгоритме триангуляции матрицы, рассмотренном ранее, можно распараллелить два цикла (один из которых является вложенным) и функцию col_max, выполняющую поиск максимума в столбце матрицы:
    . Мы будем распараллеливать с OpenMP реализацию этого алгоритма на С++.

    Для тестирования производительности параллельной программы будем использовать функцию триангуляцию следующим образом:

    #include "matrix.h"
    #include <iostream>
    #include <vector>
    #include "omp.h"
    
    int main() {
      const int n = 1000;
      std::vector<std::vector<double> > matrix(n);
    
      for (int i = 0; i < n; ++i) {
        matrix[i].resize(n);
        for (int j = 0; j < n; ++j)
          matrix[i][j] = rand();
      }
    
      double start_time = omp_get_wtime();
    
      triangulation(matrix, n);
    
      std::cout << omp_get_wtime() - start_time;
    }

    Тут создается матрица, засекается время, затем выполняется триангуляция и выбирается еще одна отсечка времени. Разница отсечек выводится на экран.

    В этой программе можно распараллелить несколько различных участков кода:

    1. поиск максимального элемента в столбце (функция col_max);
    2. обнуление с помощью эквивалентных преобразований элементов i-того столбца матрицы;
    3. добавление к j-той строке матрицы элементов i-той строки, домноженных на константу;

    В матрице находится N*N элементов. Поиск максимального элемента столбца выполняется N раз, при этом асимптотическая сложность алгоритмаO(N). После вычисления максимума выполняется перестановка строк, если она будет производиться поэлементно – трудоемкость также составит O(N), однако при использовании семантики перемещения (ее поддерживает стандартный класс vector) сложность будет всего лишь O(1). После перестановки строк выполняется обнуление элементов i-того столбца, расположенных ниже i-той строки путем эквивалентных преобразований. При этом обрабатывается N-i строк, выполняется умножение строки на константу и сложение строк (N-i элементов в строках). Трудоемкость эквивалентных преобразований составляет O(N^2), а общая сложность алгоритма – O(N^3).

    Очевидно, распараллеливание поиска максимального элемента столбца не окажет существенного влияние на время выполнения. Тем не менее, проверим это. Между потоками можно разделить только обработку элементов столбца (разные потоки должны обрабатывать разные элементы), в результате вычислений должен быть сформирован максимальный элемент столбца и его индекс. Разные потоки не должны одновременно пытаться изменять значение максимального элемента или пытаться читать его в то время, когда другой пишет. Если бы мы поместили всю работу с максимальным элементом в критическую секцию – функция фактически осталась бы последовательной, т.к. в один времени работал бы только один поток, а остальные находились бы в состоянии ожидания. В связи с этим, создадим в каждом потоке свою локальную копию максимума и вычислим максимум в части массива, соответствующей потоку. Затем, для получения максимума всего столбца выполним сравнение локальной копии с общим максимумом в критической секции. Таким образом, например для столбца из 1000 элементов при четырех потоках, каждый поток обработает по 250 элементов, а затем “последовательно” выполнится обработка четырех локальных максимумов:

    template <typename T>
    int col_max(const std::vector<std::vector<T> > &matrix, int col, int n) {
      T max = std::abs(matrix[col][col]);
      int maxPos = col;
    #pragma omp parallel
      {
        T loc_max = max;
        T loc_max_pos = maxPos;
      #pragma omp for
        for (int i = col+1; i < n; ++i) {
          T element = std::abs(matrix[i][col]);
          if (element > loc_max) {
            loc_max = element;
            loc_max_pos = i;
          }
        }
      #pragma omp critical
        {
          if (max < loc_max) {
            max = loc_max;
            maxPos = loc_max_pos;
          }
        }
      }
      return maxPos;
    }

    Распараллеливание эквивалентных преобразований реализуется проще. Между потоками можно разделить каждую из N-i строк, для которых выполняются преобразования:

        #pragma omp parallel for
        for (int j = i + 1; j < n; ++j) {
          T mul = -matrix[j][i]/matrix[i][i];
    
          for (int k = i; k < n; ++k) {
            matrix[j][k] += matrix[i][k]*mul;
          }
        }

    Либо распараллелить обработку элементов этих строк (умножение на константу и сложение):
        for (int j = i + 1; j < n; ++j) {
          T mul = -matrix[j][i]/matrix[i][i];
          #pragma omp parallel for
          for (int k = i; k < n; ++k) {
            matrix[j][k] += matrix[i][k]*mul;
          }
        }

    С точки зрения асимптотической сложности, оба варианта эквиваленты, однако во втором теоретически (может зависеть от реализации и оптимизации компилятором) затраты на создание и уничтожение потоков гораздо больше.

    Результаты тестирования распараллеливания сведены в таблицу:

    500 1000 1500 2000
    последовательная 0.509422 3.82284 12.917 30.8544
    Параллельная (col_max) 0.489255 3.86684 13.243 30.6947
    Параллельная (triangulation,
    обработка столбца)
    0.2647 2.05472 7.35282 17.202
    Параллельная (triangulation,
    обработка строки)
    0.335515 2.55324 8.79488 17.888

    Работа программы проверялась на процессоре с двумя ядрами Intel(R) Pentium(R) CPU G2130 @ 3.20GHz и шестью гигабайтами оперативной памяти. Использовался компилятор gcc 4.8.5 с опцией оптимизации -O1. Программа запускалась в операционной системе openSUSE 42.1 (x86_64).

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