Волновые схемы параллельных вычислений

Теперь рассмотрим возможность построения параллельного алгоритма, выполняющего только те вычислительные действия, что и последовательный метод и, как результат, обеспечивал бы получение точно таких же решений исходной вычислительной задачи. В последовательном алгоритме каждое очередное k -ое приближение значения вычисляется по последнему k -ому приближению значений и и предпоследнему (k-1)-ому приближению значений и . В результате при требовании совпадения результатов вычислений последовательных и параллельных вычислительных схем в начале каждой итерации метода только одно значение может быть пересчитано (невозможно распараллелить). После пересчета вычисления могут выполняться уже в двух узлах сетки и (в этих узлах выполняются условия последовательной схемы), затем после пересчета узлов и - в узлах , и и т.д. Выполнение итерации метода сеток можно разбить на последовательность шагов, на каждом из которых к вычислениям окажутся подготовленными узлы вспомогательной диагонали сетки с номером, определяемом номером этапа(рис. 6.8). Такая схема названа волной или фронтом вычислений, а алгоритмы, получаемые на ее основе, - методами волновой обработки данных (wavefront or hyperplane methods).

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

 

Таблица 6.3. Результаты экспериментов для параллельных вариантов алгоритма Гаусса-Зейделя с волновой схемой расчета (p =4)

(k – количество итераций, t – время в сек., S – ускорение)

 

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

// Алгоритм 6.5 omp_lock_t dmax_lock; omp_init_lock(dmax_lock); do { dmax = 0; // максимальное изменение значений u // нарастание волны (nx – размер волны) for (nx=1; nx<N+1; nx++) { dm[nx] = 0; #pragma omp parallel for shared(u,nx,dm) private(i,j,temp,d) for (i=1; i<nx+1; i++) {    j = nx + 1 – i;    temp = u[i][j];   u[i][j] = 0.25*(u[i-1][j]+u[i+1][j]+   u[i][j-1]+u[i][j+1]–h*h*f[i][j]); d = fabs(temp-u[i][j]);   if (dm[i] < d) dm[i] = d; } // конец параллельной области } // затухание волны for (nx=N-1; nx>0; nx--) { #pragma omp parallel for shared(u,nx,dm) private(i,j,temp,d) for (i=N-nx+1; i<N+1; i++) {    j = 2*N - nx – I + 1;    temp = u[i][j];   u[i][j] = 0.25*(u[i-1][j]+u[i+1][j]+   u[i][j-1]+u[i][j+1]–h*h*f[i][j]);   d = fabs(temp-u[i][j])   if (dm[i] < d) dm[i] = d; } // конец параллельной области } #pragma omp parallel for shared(n,dm,dmax) private(i) for (i=1; i<nx+1; i++) { omp_set_lock(dmax_lock); if (dmax < dm[i]) dmax = dm[i]; omp_unset_lock(dmax_lock); } // конец параллельной области } while (dmax > eps);

 

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

После обработки всех элементов волны среди массива dm находится максимальная погрешность выполненной итерации вычислений. Именно эта последняя часть расчетов может оказаться наиболее неэффективной из-за высоких дополнительных затрат на синхронизацию. Улучшение ситуации может быть достигнуто за счет увеличения размера последовательных участков и сокращения количества необходимых взаимодействий параллельных участков вычислений. Возможный вариант реализации такого подхода:

chunk = 200; // размер последовательного участка #pragma omp parallel for shared(n,dm,dmax) private(i,d) for (i=1; i<nx+1; i+=chunk) { d = 0; for (j=i; j<i+chunk; j++) if (d < dm[j]) d = dm[j]; omp_set_lock(dmax_lock); if (dmax < d) dmax = d; omp_unset_lock(dmax_lock); } // конец параллельной области

 

Прием укрупнения последовательных участков вычислений для снижения затрат на синхронизацию называют фрагментированием (chunking). Результаты экспериментов для данного варианта параллельных вычислений приведены в табл. 6.3.

Следует обратить внимание на еще один момент при анализе эффективности разработанного параллельного алгоритма. Фронт волны вычислений плохо соответствует правилам использования кэша - быстродействующей дополнительной памяти компьютера, используемой для хранения копии наиболее часто используемых областей оперативной памяти. Использование кэша может существенно повысить (в десятки раз) быстродействие вычислений. Размещение данных в кэше может происходить предварительно (при использовании алгоритмов предсказания потребности в данных) или в момент извлечения значений из основной оперативной памяти. Подкачка данных в кэш осуществляется не одиночными значениями, а небольшими группами – строками кэша (cache line). Загрузка значений в строку кэша осуществляется из последовательных элементов памяти; типовые размеры строки кэша обычно равны 32, 64, 128, 256 байтам. Эффект наличия кэша будет виден, если выполняемые вычисления используют одни и те же данные многократно (локальность обработки данных) и осуществляют доступ к элементам памяти с последовательно возрастающими адресами (последовательность доступа).

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

Рис. 6.9. Блочное представление сетки области расчетов

 

Создаваемый на основе такого подхода метод вычислений в общем виде может быть описан следующим образом (блоки образуют в области расчётов прямоугольную решётку размера NBxNB):

// Алгоритм 6.6 // NB количество блоков do { // нарастание волны (размер волны равен nx+1) for (nx =0; nx<NB; nx++) { #pragma omp parallel for shared(nx) private(i,j) for (i=0; i<nx+1; i++) {    j = nx – i;    // <обработка блока с координатами (i,j)> } // конец параллельной области } // затухание волны for (nx=NB-2; nx>-1; nx--) { #pragma omp parallel for shared(nx) private(i,j) for (i=0; i<nx+1; i++) {    j = 2*(NB-1) - nx – i;    // <обработка блока с координатами (i,j)> } // конец параллельной области } // <определение погрешности вычислений> } while (dmax > eps);

 

Вычисления в предлагаемом алгоритме происходят в соответствии с волновой схемой обработки данных – вначале вычисления выполняются только в левом верхнем блоке с координатами (0,0), далее для обработки становятся доступными блоки с координатами (0,1) и (1,0) и т.д. – см. результаты экспериментов в табл. 6.3.

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

Наилучшие показатели использования кэша будут достигаться, если в кэше будет достаточно места для размещения не менее трех строк блока (при обработке строки блоки блока используются данные трех строк блока одновременно). Тем самым, исходя из размера кэша, можно определить рекомендуемый максимально-возможный размер блока. Так, например, при кэше 8 Кб и 8-байтовых значениях данных этот размер составит приближенно 300 (8Кб/3/8). Можно определить и минимально-допустимый размер блока из условия совпадения размеров строк кэша и блока. Так, при размере строки кэша 256 байт и 8-байтовых значениях данных размер блока должен быть кратен 32.

Последнее замечание следует сделать о взаимодействии граничных узлов блоков. Учитывая граничное взаимодействие, соседние блоки целесообразно обрабатывать на одних и тех же процессорах. В противном случае, можно попытаться так определить размеры блоков, чтобы объем пересылаемых между процессорами граничных данных был минимален. Так, при размере строки кэша в 256 байт, 8-байтовых значениях данных и размере блока 64х64 объем пересылаемых данных 132 строки кэша, при размере блока 128х32 – всего 72 строки. Такая оптимизация имеет наиболее принципиальное значение при медленных операциях пересылки данных между кэшами процессоров, т.е. для систем с неоднородным доступом к памяти, (nonuniform memory access - NUMA).