Эффективное распараллеливание линейной алгебраической функции в C++ OpenMP
У меня мало опыта в параллельном программировании, и мне было интересно, если бы кто-нибудь мог бы бегло взглянуть на немного кода, который я написал, и посмотреть, есть ли какие-нибудь очевидные способы, которыми я могу повысить эффективность вычислений.
Трудность возникает из-за того, что у меня есть несколько матричных операций неравномерной размерности, которые мне нужно вычислить, поэтому я не уверен, что это наиболее сжатый способ кодирования вычислений.
Ниже мой код. Обратите внимание, этот код работает. Матрицы, с которыми я работаю, имеют размер приблизительно 700x700 (см. Ниже) или 700x30 [int n].
Кроме того, я использую библиотеку броненосца для моего последовательного кода. Может случиться так, что распараллеливание с использованием openMP, но с сохранением классов матрицы броненосца происходит медленнее, чем по умолчанию для стандартной библиотеки; У кого-нибудь есть мнение по этому поводу (до того, как я потрачу часы на капитальный ремонт!)?
double start, end, dif;
int i,j,k; // iteration counters
int s,n; // matrix dimensions
mat B; B.load(...location of stored s*n matrix...) // input objects loaded from file
mat I; I.load(...s*s matrix...);
mat R; R.load(...s*n matrix...);
mat D; D.load(...n*n matrix...);
double e = 0.1; // scalar parameter
s = B.n_rows; n = B.n_cols;
mat dBdt; dBdt.zeros(s,n); // object for storing output of function
// 100X sequential computation using Armadillo linear algebraic functionality
start = omp_get_wtime();
for (int r=0; r<100; r++) {
dBdt = B % (R - (I * B)) + (B * D) - (B * e);
}
end = omp_get_wtime();
dif = end - strt;
cout << "Seq computation: " << dBdt(0,0) << endl;
printf("relaxation time = %f", dif);
cout << endl;
// 100 * parallel computation using OpenMP
omp_set_num_threads(8);
for (int r=0; r<100; r++) {
// parallel computation of I * B
#pragma omp parallel for default(none) shared(dBdt, B, I, R, D, e, s, n) private(i, j, k) schedule(static)
for (i = 0; i < s; i++) {
for (j = 0; j < n; j++) {
for (k = 0; k < s; k++) {
dBdt(i, j) += I(i, k) * B(k, j);
}
}
}
// parallel computation of B % (R - (I * B))
#pragma omp parallel for default(none) shared(dBdt, B, I, R, D, e, s, n) private(i, j) schedule(static)
for (i = 0; i < s; i++) {
for (j = 0; j < n; j++) {
dBdt(i, j) = R(i, j) - dBdt(i, j);
dBdt(i, j) *= B(i, j);
dBdt(i, j) -= B(i, j) * e;
}
}
// parallel computation of B * D
#pragma omp parallel for default(none) shared(dBdt, B, I, R, D, e, s, n) private(i, j, k) schedule(static)
for (i = 0; i < s; i++) {
for (j = 0; j < n; j++) {
for (k = 0; k < n; k++) {
dBdt(i, j) += B(i, k) * D(k, j);
}
}
}
}
end = omp_get_wtime();
dif = end - strt;
cout << "OMP computation: " << dBdt(0,0) << endl;
printf("relaxation time = %f", dif);
cout << endl;
Если я использую Hyper-Threading 4 ядра, я получаю следующий вывод:
Seq computation: 5.54926e-10
relaxation time = 0.130031
OMP computation: 5.54926e-10
relaxation time = 2.611040
Это говорит о том, что хотя оба метода дают один и тот же результат, параллельная формулировка примерно в 20 раз медленнее, чем последовательная.
Возможно, что для матриц такого размера издержки, связанные с этой проблемой "переменной размерности", перевешивают преимущества распараллеливания. Любое понимание будет высоко ценится.
Заранее спасибо,
Джек
2 ответа
Искусство HPC- это правильная эффективность (плохие гранты никогда не получают кластерную квоту HPC)
- так что первая надежда - ваш процесс никогда не будет перечитан из файла
Зачем? Это будет HPC-убийца:
Мне нужно повторить это вычисление много тысяч раз
Достаточно справедливо сказать, что этот комментарий увеличил общую потребность в полном пересмотре подхода и переделке будущего решения не для того, чтобы полагаться на несколько уловок, а для того, чтобы действительно извлечь выгоду из вашего конкретного случая.
Не в последнюю очередь - [PARALLEL] планирование не нужно, так как "просто" - [CONCURRENT]-Процесс планирования здесь вполне достаточно. Нет необходимости организовывать какую-либо явную межпроцессную синхронизацию или передачу сообщений, и процесс может быть просто организован для достижения наилучшей возможной производительности.
Никакой "... быстрый взгляд на немного кода..." не поможет
Сначала вы должны понять как весь процесс, так и аппаратные ресурсы, на которых он будет выполняться.
Тип процессора покажет вам доступные расширения набора команд для продвинутых трюков, размеры кеша L3- / L2- / L1 + размеры строк кеша помогут вам выбрать наилучшее дружественное к кешу повторное использование дешевого доступа к данным (не платя) сотни [нс], если вместо этого можно работать умнее всего на несколько [нс] из еще не выселенной NUMA-core-local копии)
Математика первая, реализация следующая:
Как дано dBdt = B % ( R - (I * B) ) + ( B * D ) - ( B * e )
При ближайшем рассмотрении каждый должен быть готов реализовать приоритеты HPC/ выравнивания кэша и ошибочно замкнутые ловушки:
dBdt = B % ( R - ( I * B ) ) ELEMENT-WISE OP B[s,n]-COLUMN-WISE
+ ( B * D ) SUM.PRODUCT OP B[s,n].ROW-WISE MUL-BY-D[n,n].COL
- ( B * e ) ELEMENT-WISE OP B[s,n].ROW-WISE MUL-BY-SCALAR
ROW/COL-SUM.PRODUCT OP -----------------------------------------+++++++++++++++++++++++++++++++++++++++++++++
ELEMENT-WISE OP ---------------------------------------------+ |||||||||||||||||||||||||||||||||||||||||||||
ELEMENT-WISE OP ----------------------+ | |||||||||||||||||||||||||||||||||||||||||||||
| | |||||||||||||||||||||||||||||||||||||||||||||
v v vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
dBdt[s,n] = B[s,n] % / R[s,n] - / I[s,s] . B[s,n] \ \
_________[n] _________[n] | _________[n] | ________________[s] _________[n] | |
|_| | |_| | | |_| | | |________________| | | | | |
| . | | . | | | . | | | | | | | | |
| . | | . | | | . | | | | | | | | |
| . | | . | | | . | | | | | | | | |
| . | = | . | % | | . | - | | | . | | | | |
| . | | . | | | . | | | | | | | | |
| . | | . | | | . | | | | | | | | |
| . | | . | | | . | | | | | | | | |
[s]|_________| [s]|_________| | [s]|_________| | [s]|________________| [s]|_|_______| | |
\ \ / /
B[s,n] D[n,n]
_________[n] _________[n]
|_________| | | |
| . | | | |
| . | | | |
| . | | | |
+ | . | . [n]|_|_______|
| . |
| . |
| . |
[s]|_________|
B[s,n]
_________[n]
|_| . . . |
| . |
| . |
| . |
- | . | * REGISTER_e
| . |
| . |
| . |
[s]|_________|
Имея это в виду, эффективные петли HPC будут сильно отличаться
В зависимости от реальных CPU-кешей,
цикл может очень эффективно совместно обрабатывать B -строка выровнен ( B * D ) - ( B * e )
в одной фазе, а также на основе самой высокой эффективности повторного использования поэлементно самого длинного трубопровода B % ( R - ( I * B ) ) здесь есть возможность повторно использовать ~ 1000 x ( n - 1) попаданий в кэш B - выровненный по столбцам, который вполне должен вписаться в следы кеша L1-DATA, таким образом, достигается экономия порядка нескольких секунд только благодаря циклам, выровненным по кешу.
После того, как это дружественное к кешу выравнивание цикла закончено,
Далее может помочь распределенная обработка, а не раньше.
Итак, план эксперимента:
Шаг 0: Основа-Истина: ~ 0.13 [s] за dBdt[700,30] используя броненосец в 100-тестовых петлях
Шаг 1. Последовательность руководства: - проверить награды за лучший код с выравниванием кеша (не отправленный, а математически эквивалентный, оптимизированный для повторного использования строк кэша), где должно быть не более 4х for(){...} 2-вложенные кодовые блоки, остальные 2 внутри, чтобы соответствовать правилам линейной алгебры, без разрушительных преимуществ выравнивания строк кеша (с некоторым остаточным потенциалом, чтобы получить выгоду еще больше в [PTIME] от использования дублированного [PSPACE]расположение данных (как в FORTRAN-порядке, так и в C-порядке для соответствующих стратегий повторного чтения), поскольку матрицы имеют настолько миниатюрные размеры, а L2- / L1-DATA-кэш-память, доступная для каждого ядра ЦП, имеет размеры кеша, хорошо растущие в масштаб)
Шаг 2: Руководство omp (<= NUMA_cores - 1): - проверить, действительно ли omp может привести к какому-либо "положительному" ускорению закона Амдала (помимо накладных расходов на установку omp). Тщательное сопоставление аффинности process-2-CPU_core может помочь избежать любого возможного вытеснения кеша, созданного любым не-HPC-потоком, портящим дружественный кеш-макет на наборе "зарезервированных" наборов конфигурации ( NUMA_cores - 1 ) где все остальные (не-HPC-процессы) должны отображаться на последнем (совместно используемом) ядре ЦП, что помогает предотвратить сохранение этими строками-кеш-ядрами своих строк кэша, не удаленных ни одним ядром / планировщиком. не-HPC-нить.
(Как видно из (2), существуют рекомендации, основанные на передовых практиках HPC, которые ни один компилятор (даже оснащенный волшебной палочкой) никогда не сможет реализовать, поэтому не стесняйтесь обратиться к своему PhD Tutor за помощью. с другой стороны, если ваш тезис нуждается в некотором HPC-опыте, поскольку в этой довольно дорогой экспериментальной области не так-то просто построить пробную ошибку, а вашим основным доменом является не линейная алгебра и / или окончательный CS-теоретический / HW-специфичный кеш -стратегия оптимизаций.)
Эпилог:
Использование интеллектуальных инструментов ненадлежащим образом не приводит к чему-то большему, чем дополнительные накладные расходы (разбиение / объединение задач + трансляция памяти (хуже с атомарной блокировкой (хуже с блокировкой / ограждением / барьерами))).
Если вы используете компилятор, который исправляет гнезда плохих циклов и циклы плавких предохранителей, чтобы улучшить локальность памяти для непараллельных сборок, openmp, вероятно, отключит эту оптимизацию. В соответствии с рекомендациями других, вы должны рассмотреть оптимизированную библиотеку, такую как mkl или acml. Стандартные gfortran blas, обычно поставляемые с дистрибутивами, не являются многопоточными.