Русский
Русский
English
Статистика
Реклама

Немного об ускорении программы распараллеливание (ручное или автоматическое) на базе сверхоптимистичных вычислений

Здравствуйте, уважаемые читатели. В этой публикации речь пойдет о такой (уже ставшей привычной) вещи как ускорение работы программы путем применения параллельных вычислений. Технологии организации таких вычислений известны это и обычное многопоточное программирование, и применение специальных интерфейсов: OpenMP, OpenAcc, MPI, DVM и многих других (при этом распараллеливаются циклы, используется векторизация или конвейеризация, организуются ленивые вычисления, выделяются независимые блоки программы, которые можно запустить в параллель и т.п.).

При этом обычно исходят из той идеи, что распараллеливание не должно каким-то образом влиять на результаты исполнения программы. Это жесткое, но справедливое для многих случаев требование. Однако если мы пытаемся распараллелить программу, ведущую какие-либо расчеты численными методами (обучаем нейронную сеть, моделируем динамику жидкости или молекулярной системы, решаем обыкновенные дифференциальные уравнения или оптимизационные задачи), то результат и так (в любом случае) будет иметь некоторую погрешность. Поэтому, почему бы не применить рискованные технологии распараллеливания, которые могут внести в математическое решение небольшую дополнительную погрешность, но позволят получить еще некоторое дополнительное ускорение? Об одной из таких технологий о расщеплении тел циклов с предсказанием промежуточных результатов и откатом при неудачном предсказании (собственно, это и есть сверхоптимистичные вычисления в частично транзакционной памяти) и пойдет речь.

Идея распараллеливания


Предположим, что мы имеем цикл, тело которого состоит из двух последовательных частей, причем вторая часть зависит от первой. Пусть и отдельные витки цикла зависят друг от друга. Например:

for (int i = 0; i < N; i++) {x = f(y);y = g(x);}

На первый взгляд, распараллелить такой цикл невозможно. Однако мы попробуем. Попытаемся исполнять параллельно первый и второй операторы тела цикла. Проблема состоит в том, что на момент вычисления g(x) должен быть известен x, но он будет рассчитан только в конце первой части. Что же, введем некоторую схему, которая в начале второй части попытается предсказать новое значение x. Можно это сделать, например, с помощью линейной предикции, которая обучится предсказывать новое значение x, опираясь на историю его изменения. Тогда вторую часть можно считать параллельно с первой (это и есть сверхоптимизм), а когда обе будут подсчитаны, сравнить предсказанное значение x с реальным, полученным в конце первой части. Если они примерно равны, то результат вычислений обеих частей можно принять (и перейти к следующему витку цикла). А если они сильно отличаются, то потребуется пересчитать только вторую часть. При такой схеме в какой-то части случаев получим чистое распараллеливание, в остальных фактический последовательный счет. Алгоритм выполнения цикла при этом такой:

for (int i = 0; i < N; i++) {Распараллеливаем на два ядра {На ядре 1  считаем x = f(y). Далее передаем во вторую часть получение значение x;На ядре 2  предсказываем значение x* и считаем y* = g(x*). Получаем значение x из первой части и сравниваем его с x*. Если разница невелика, то y = y* и завершаем итерацию цикла. Если различие большое, повторяем вычисление с новыми данными: y = g(x). }}

Базовый алгоритм ясен. Теоретическое ускорение в два раза, но на практике будет, конечно, меньше, поскольку: а) часть времени тратится на предсказания и согласования; б) не все итерации выполнятся параллельно; в) первая и вторая части тела цикла могут иметь различную трудоемкость (в идеале требуется равная). Перейдем к реализации.

Реализация распараллеливания сверхоптимистичные вычисления


Поскольку в алгоритме распараллеливания идет речь об отмене части расчетов (при неудаче) и их повторном выполнении, здесь явно есть что-то от идеи работы в транзакционной памяти. Лучше в частично транзакционной, где определенные переменные работают по схеме транзакционной памяти, а остальные переменные как обычно. Передачу данных из первой части во вторую можно организовать с помощью некоторого специального канала. Пусть этот канал будет предсказывающим: а) если на момент приема данные в канал уже переданы, то они из него и читаются, б) если на момент приема данные в канал еще не поступили, то он пытается предсказать эти данные и возвращает результат предсказания. Этот канал и будет работать по немного не свойственной обычной транзакционной памяти схеме: если в конце транзакции второй части цикла обнаружится расхождение между поступившими в канал данными и предсказанными им данными, то транзакция второй части цикла отменяется и исполняется повторно, при этом из канала будут читаться уже не предсказания, а реально пришедшие данные. Цикл приобретет вид:

for (int i = 0; i < N; i++) {Распараллеливаем на два ядра, включаем частично транзакционную память {Ядро 1 (транзакция 1):x = f(y);Предсказывающий_Канал.put(x);Ядро 2 (транзакция 2):Предсказывающий_Канал.get(x);y = g(x);}}

Готово. Заботу о предсказании данных взял на себя канал, заботу об отмене расчетов при излишне оптимистичном предсказании взяла на себя частично транзакционная память.

Некоторые полезные применения: нейронные сети, метод частиц в ячейках


Такую схему распараллеливания тела цикла можно применить в ряде математических алгоритмов, например, при моделировании электростатической линзы методом частиц в ячейках, а также при обучении нейронной сети прямого распространения методом обратного распространения ошибки. Первая задача очень специальная, поэтому обсуждать ее здесь я не буду, скажу только, что изложенный подход к распараллеливанию дал ускорение на 10-15%. А вот вторая задача уже более популярная, поэтому о ней несколько слов сказать просто необходимо.

Цикл обучения нейронной сети включает последовательный проход по обучающим парам, причем для каждой пары выполняется прямой ход (расчет выхода сети) и обратный ход (коррекция весов и смещений). Это и есть две части тела цикла по обучающим парам и для их распараллеливания можно применить вышеизложенный подход (кстати, его можно применить и при параллельном проходе по обучающим парам, с незначительными изменениями). В результате на типичной задаче обучения нейронной сети я получил 50% выигрыша по скорости работы.

Автоматизация распараллеливания C-программ


Идея сверхоптимистичных вычислений не очень сложна, поэтому была написана специальная программа-транслятор, которая занимается автоматическим распараллеливанием находит в исходной C-программе циклы, для которых такое распараллеливание может дать положительный результат и расщепляет их тела на две части, вставляя необходимые директивы OpenMP, находя потенциальные переменные для каналов, подключая библиотеку работы с частично транзакционной памятью и предицирующими каналами и, в конечном итоге, порождая выходную распараллеленную программу.

В частности, такой транслятор был применен к программе моделирования электростатической линзы. Приведу обе программы исходную (в которую включена директива-указание на распараллеливание циклов) и полученную после трансляции.

Исходная программа:


#include <stdlib.h>#include <stdio.h>#include <math.h>#pragma auto parallelize#pragma auto pure(malloc,fabs,free,sizeof,omp_get_wtime)#define theta 1.83#define NX 40#define NY 40#define h 0.1#define NP 15000// Собирающая электростатическая линза#define U1 200#define U2 5000#define e -1.5E-13#define m 1E-11#define e0 8.85E-12#define V (h*h)#define tau 0.000015#define T 0.09#define POISSON_EPS 0.01#define TOL_EPS 0.25int main() {        double * U  = (double *)malloc(NY*NX*sizeof(double));        double * UU = (double *)malloc(NY*NX*sizeof(double));        double * EX = (double *)malloc(NY*NX*sizeof(double));        double * EY = (double *)malloc(NY*NX*sizeof(double));double * PX = (double *)malloc(NP*sizeof(double));double * PY = (double *)malloc(NP*sizeof(double));int * X = (int *)malloc(NP*sizeof(int));int * Y = (int *)malloc(NP*sizeof(int));double ro[NY][NX];split_private double t;split_private double tm;split_private int i, j;for (i = 0; i < NY; i++)for (j = 0; j < NX; j++) {UU[i*NX+j] = j == NX-1 ? U2 : j == NX/2 && (i < NY/4 || i > 3*NY/4) ? U1 : 0.0;EX[i*NX+j] = 0.0;EY[i*NX+j] = 0.0;}for (i = 0; i < NP; i++) {int x, y;PX[i] = 0.5*NX*h*rand()/RAND_MAX;PY[i] = NY*h*rand()/RAND_MAX;x = PX[i]/h;y = PY[i]/h;if (x < 0) x = 0;else if (x > NX-1) x = NX-1;if (y < 0) y = 0;else if (y > NY-1) y = NY-1;X[i] = x;Y[i] = y;}tm = omp_get_wtime();for (t = 0.0; t < T; t += tau) {unsigned int n[NY][NX] = { 0 };double err;int ptr = 0;for (i = 0; i < NY; i++)    for (j = 0; j < NX; j++, ptr++)U[ptr] = UU[ptr];for (i = 1; i < NY - 1; i++)for (j = 1; j < NX - 1; j++) {EX[i*NX+j] = -(U[i*NX+j+1]-U[i*NX+j-1])/2.0/h;EY[i*NX+j] = -(U[(i+1)*NX+j]-U[(i-1)*NX+j])/2.0/h;}for (i = 0; i < NP; i++) {PX[i] += tau*e*EX[Y[i]*NX+X[i]]/m;PY[i] += tau*e*EY[Y[i]*NX+X[i]]/m;}for (i = 0; i < NP; i++) {int x = PX[i]/h;int y = PY[i]/h;if (x < 0) x = 0;else if (x > NX-1) x = NX-1;if (y < 0) y = 0;else if (y > NY-1) y = NY-1;Y[i] = y;X[i] = x;n[y][x]++;}for (i = 0; i < NY; i++)for (j = 0; j < NX; j++)ro[i][j] = n[i][j]*e/V;do {err = 0.0;for (i = 1; i < NY - 1; i++)for (j = 1+(i-1)%2; j < NX - 1; j+=2) {  int ptr = i*NX + j;  if (!(j == NX/2 && (i < NY/4 || i > 3*NY/4))) {double _new = (1-theta)*UU[ptr] + theta/4.0*(UU[ptr-1]+UU[ptr+1]+UU[ptr+NX]+UU[ptr-NX]-h*h*ro[i][j]/e0);double loc_err = fabs(UU[ptr] - _new);if (loc_err > err) err = loc_err;UU[ptr] = _new;  }}for (i = 1; i < NY - 1; i++)for (j = 1+i%2; j < NX - 1; j+=2) {  int ptr = i*NX + j;  if (!(j == NX/2 && (i < NY/4 || i > 3*NY/4))) {double _new = (1-theta)*UU[ptr] + theta/4.0*(UU[ptr-1]+UU[ptr+1]+UU[ptr+NX]+UU[ptr-NX]-h*h*ro[i][j]/e0);double loc_err = fabs(UU[ptr] - _new);if (loc_err > err) err = loc_err;UU[ptr] = _new;  }}for (j = 0; j < NX; j++) {UU[j] = UU[NX + j];UU[(NY-1)*NX + j] = UU[(NY-2)*NX + j];}} while (err > POISSON_EPS);}for (i = 0; i < NY; i++) {for (j = 0; j < NX; j++)printf("%lf\t", UU[i*NX+j]);printf("\n");}return 0;}

Автоматически распараллеленная программа


#include "transact.h"#define split_private /* split-private */#include <stdlib.h>#include <stdio.h>#include <math.h>#define theta 1.83#define NX 40#define NY 40#define h 0.1#define NP 15000#define U1 200#define U2 5000#define e -1.5E-13#define m 1E-11#define e0 8.85E-12#define V (h*h)#define tau 0.000015#define T 0.09#define POISSON_EPS 0.01#define TOL_EPS 0.25int  main(  ){  double * U  = (double *)malloc(NY*NX*sizeof(double));  double * UU = (double *)malloc(NY*NX*sizeof(double));  double * EX = (double *)malloc(NY*NX*sizeof(double));  double * EY = (double *)malloc(NY*NX*sizeof(double));  double * PX = (double *)malloc(NP*sizeof(double));  double * PY = (double *)malloc(NP*sizeof(double));  int * X = (int *)malloc(NP*sizeof(int));  int * Y = (int *)malloc(NP*sizeof(int));  double ro[NY][NX];  split_private double t;  split_private double tm;  split_private int i, j;  for ( i = 0; i < NY; i++ )    for ( j = 0; j < NX; j++ )      {        UU[i*NX+j] = j == NX-1 ? U2 : j == NX/2 && (i < NY/4 || i > 3*NY/4) ? U1 : 0.0;        EX[i*NX+j] = 0.0;        EY[i*NX+j] = 0.0;      }  for ( i = 0; i < NP; i++ )    {      int x, y;      PX[i] = 0.5*NX*h*rand()/RAND_MAX;      PY[i] = NY*h*rand()/RAND_MAX;      x = PX[i]/h;      y = PY[i]/h;      if ( x < 0 )        x = 0;      else        if ( x > NX-1 )          x = NX-1;      if ( y < 0 )        y = 0;      else        if ( y > NY-1 )          y = NY-1;      X[i] = x;      Y[i] = y;    }  tm = omp_get_wtime();#pragma omp parallel num_threads(2) private(t,tm,i,j)   {    int __id__ = omp_get_thread_num();    TOut<double > * out_ro = __id__ == 0 ? new TOut<double >("ro63", (NY)*(NX), 2, 0.01, -1, "63") : NULL;    TIn<double > * in_ro = __id__ == 1 ? new TIn<double >("ro63", (NY)*(NX), 2, 0.01, -1, "63") : NULL;    for ( t = 0.0; t < T; t += tau )      {        unsigned int n[NY][NX] = { 0 };        double err;        int ptr = 0;        if ( __id__ == 0 )          {            for ( i = 0; i < NY; i++ )              for ( j = 0; j < NX; j++, ptr++ )                U[ptr] = UU[ptr];          }transaction_atomic("63")        {          if ( __id__ == 0 )            {              for ( i = 1; i < NY - 1; i++ )                for ( j = 1; j < NX - 1; j++ )                  {                    EX[i*NX+j] = -(U[i*NX+j+1]-U[i*NX+j-1])/2.0/h;                    EY[i*NX+j] = -(U[(i+1)*NX+j]-U[(i-1)*NX+j])/2.0/h;                  }              for ( i = 0; i < NP; i++ )                {                  PX[i] += tau*e*EX[Y[i]*NX+X[i]]/m;                  PY[i] += tau*e*EY[Y[i]*NX+X[i]]/m;                }              for ( i = 0; i < NP; i++ )                {                  int x = PX[i]/h;                  int y = PY[i]/h;                  if ( x < 0 )                    x = 0;                  else                    if ( x > NX-1 )                      x = NX-1;                  if ( y < 0 )                    y = 0;                  else                    if ( y > NY-1 )                      y = NY-1;                  Y[i] = y;                  X[i] = x;                  n[y][x]++;                }              for ( i = 0; i < NY; i++ )                for ( j = 0; j < NX; j++ )                  ro[i][j] = n[i][j]*e/V;              out_ro->put((double  *)ro);            }          else            {              double  ro[NY][NX];              in_ro->get((double  *)ro, 0);              do                {                  err = 0.0;                  for ( i = 1; i < NY - 1; i++ )                    for ( j = 1+(i-1)%2; j < NX - 1; j+=2 )                      {                        int ptr = i*NX + j;                        if ( !(j == NX/2 && (i < NY/4 || i > 3*NY/4)) )                          {                            double _new = (1-theta)*UU[ptr] + theta/4.0*(UU[ptr-1]+UU[ptr+1]+UU[ptr+NX]+UU[ptr-NX]-h*h*ro[i][j]/e0);                            double loc_err = fabs(UU[ptr] - _new);                            if ( loc_err > err )                              err = loc_err;                            UU[ptr] = _new;                          }                      }                  for ( i = 1; i < NY - 1; i++ )                    for ( j = 1+i%2; j < NX - 1; j+=2 )                      {                        int ptr = i*NX + j;                        if ( !(j == NX/2 && (i < NY/4 || i > 3*NY/4)) )                          {                            double _new = (1-theta)*UU[ptr] + theta/4.0*(UU[ptr-1]+UU[ptr+1]+UU[ptr+NX]+UU[ptr-NX]-h*h*ro[i][j]/e0);                            double loc_err = fabs(UU[ptr] - _new);                            if ( loc_err > err )                              err = loc_err;                            UU[ptr] = _new;                          }                      }                  for ( j = 0; j < NX; j++ )                    {                      UU[j] = UU[NX + j];                      UU[(NY-1)*NX + j] = UU[(NY-2)*NX + j];                    }                }              while ( err > POISSON_EPS )                ;            }        }      }    delete in_ro;    delete out_ro;  }  for ( i = 0; i < NY; i++ )    {      for ( j = 0; j < NX; j++ )        printf("%lf\t", UU[i*NX+j]);      printf("\n");    }  return 0;}

Итоги


Итак, иногда можно пытаться распараллелить программу даже в случаях, когда она состоит из строго последовательных фрагментов, и даже получать положительные результаты по ускорению (в моих экспериментах прирост ускорения от 15 до 50%). Надеюсь, что эта небольшая статья окажется кому-нибудь полезной.
Источник: habr.com
К списку статей
Опубликовано: 18.08.2020 20:07:52
0

Сейчас читают

Комментариев (0)
Имя
Электронная почта

Программирование

C++

Параллельное программирование

Распараллеливание

Транзакционная память

Предсказание

Язык программирования

Сверхоптимистичные вычисления

Категории

Последние комментарии

  • Имя: Макс
    24.08.2022 | 11:28
    Я разраб в IT компании, работаю на арбитражную команду. Мы работаем с приламы и сайтами, при работе замечаются постоянные баны и лаги. Пацаны посоветовали сервис по анализу исходного кода,https://app Подробнее..
  • Имя: 9055410337
    20.08.2022 | 17:41
    поможем пишите в телеграм Подробнее..
  • Имя: sabbat
    17.08.2022 | 20:42
    Охренеть.. это просто шикарная статья, феноменально круто. Большое спасибо за разбор! Надеюсь как-нибудь с тобой связаться для обсуждений чего-либо) Подробнее..
  • Имя: Мария
    09.08.2022 | 14:44
    Добрый день. Если обладаете такой информацией, то подскажите, пожалуйста, где можно найти много-много материала по Yggdrasil и его уязвимостях для написания диплома? Благодарю. Подробнее..
© 2006-2023, personeltest.ru