Декомпозиция LU с использованием openmp
У меня проблема: параллельная версия алгоритма декомпозиции LU работает в то же время, что и последовательность:
void lup_od_omp(double* a, int n){
int i,j,k;
for(k = 0; k < n - 1; ++k)
{
#pragma omp parallel for shared(a,n,k) private(i,j)
for(i = k + 1; i < n; i++)
{
a[i*n + k] /= a[k*n + k];
for(j = k + 1; j < n; j++)
{
a[i*n + j] -= a[i*n + k]*a[k*n + j];
}
}
}}
Может я что-то не так делаю?
Решение
Поскольку вы работаете только с двумя ядрами, распараллеливание может помешать векторизатору. Векторизация на SSE2 даст вам полосу пропускания данных 2 раза за операцию, 4 на AVX.
Двойной поток имеет много накладных расходов на синхронизацию, что может замедлить вас, особенно если вы потеряете векторизацию. Также по какой-то причине ваш #pragma omp не запускаются никакие темы, если omp_set_num_threads был вызван, чтобы фактически заставить его использовать потоки.
Еще одна вещь, которая также связана с векторизацией, заключается в том, что не все компиляторы понимают, что a[i*n + j] предназначен для обращения к двумерному массиву, поэтому лучше сначала объявить его таковым.
Вот небольшая оптимизация вашего кода, который довольно хорошо работает на моем Xeon:
void lup_od_omp(int n, double (*a)[n]){
int i,k;
for(k = 0; k < n - 1; ++k) {
// for the vectoriser
for(i = k + 1; i < n; i++) {
a[i][k] /= a[k][k];
}
#pragma omp parallel for shared(a,n,k) private(i) schedule(static, 64)
for(i = k + 1; i < n; i++) {
int j;
const double aik = a[i][k]; // some compilers will do this automatically
for(j = k + 1; j < n; j++) {
a[i][j] -= aik * a[k][j];
}
}
}
}
Время выполнения для массива 3000×3000 icc -O2:
Your code sequential: 0:24.61 99% CPU
Your code 8 threads : 0:05.21 753% CPU
My code sequential: 0:18.53 99% CPU
My code 8 threads : 0:05.42 766% CPU
И на другом компьютере я тестировал его на AVX (256-битные векторы, 4 двойных за операцию):
My code on AVX sequential : 0:09.45 99% CPU
My code on AVX 8 threads : 0:03.92 766% CPU
Как вы можете видеть, я немного улучшил векторизатор, но мало что сделал для параллельной секции.
Другие решения
Основная проблема вашего кода заключается в том, что вы неправильно разложили нагрузку.
Для разложения LU вы вызываете параллель для n-1 раз. В каждый раз, параллельный for будет выполнять ветвление и соединение потоков, что приводит к большим накладным расходам. Особенно когда k большой, внутренний цикл (for(i){for(j){...}}) содержит только очень мало работы. Распараллеливание будет совсем неэффективным.
Вы можете рассмотреть возможность использования надлежащих схем агломерации для сокращения накладных расходов. Для получения дополнительной информации, пожалуйста, обратитесь к этим слайдам.
http://courses.engr.illinois.edu/cs554/notes/06_lu_8up.pdf
С другой стороны, вы можете использовать существующие библиотеки производительности для получения максимальной производительности при факторизации LU, такие как Intel MKL