При этом обычно исходят из той идеи, что распараллеливание не должно каким-то образом влиять на результаты исполнения программы. Это жесткое, но справедливое для многих случаев требование. Однако если мы пытаемся распараллелить программу, ведущую какие-либо расчеты численными методами (обучаем нейронную сеть, моделируем динамику жидкости или молекулярной системы, решаем обыкновенные дифференциальные уравнения или оптимизационные задачи), то результат и так (в любом случае) будет иметь некоторую погрешность. Поэтому, почему бы не применить «рискованные» технологии распараллеливания, которые могут внести в математическое решение небольшую дополнительную погрешность, но позволят получить еще некоторое дополнительное ускорение? Об одной из таких технологий – о расщеплении тел циклов с предсказанием промежуточных результатов и откатом при неудачном предсказании (собственно, это и есть «сверхоптимистичные» вычисления в частично транзакционной памяти) и пойдет речь.
Идея распараллеливания
Предположим, что мы имеем цикл, тело которого состоит из двух последовательных частей, причем вторая часть зависит от первой. Пусть и отдельные витки цикла зависят друг от друга. Например:
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.25 int 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.25 int 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%). Надеюсь, что эта небольшая статья окажется кому-нибудь полезной.
ссылка на оригинал статьи https://habr.com/ru/post/515672/
Добавить комментарий