Простенькая задача:
Дан большой (миллиарды) массив комплексных чисел.
Надо: мин и мах по реальной и комплексной компоненте и по амплитуде.
Я на фортране решил "в лоб" с неплохой масштабируемостью (~7х на 8-ми ядрах).
Для 440 миллионов чисел: 16 секунд на одном ядре против 2.2 секунд на восьми.
Все решение заняло 15 минут на кодирование и три десятка строк.
Потом захотелось "красивостей".
Решили написать то же на С++.
Тут же обнаружили, что мин/мах редукции в С нет (ни у МS, ни у Интела).
"Изобрели" за час, но получился монстр.
Скорость примерно та же, как и в Фортране (на интеловском С).
Но выглядит преотвратно.
(К тому же, требует ОМП длл, а у Фортрана есть статическая библиотека)
Праздный вопрос (праздный, потому как слинковали с фортраном): как в С делать "красивую" редукцию?
FORTRAN
Дан большой (миллиарды) массив комплексных чисел.
Надо: мин и мах по реальной и комплексной компоненте и по амплитуде.
Я на фортране решил "в лоб" с неплохой масштабируемостью (~7х на 8-ми ядрах).
Для 440 миллионов чисел: 16 секунд на одном ядре против 2.2 секунд на восьми.
Все решение заняло 15 минут на кодирование и три десятка строк.
Потом захотелось "красивостей".
Решили написать то же на С++.
Тут же обнаружили, что мин/мах редукции в С нет (ни у МS, ни у Интела).
"Изобрели" за час, но получился монстр.
Скорость примерно та же, как и в Фортране (на интеловском С).
Но выглядит преотвратно.
(К тому же, требует ОМП длл, а у Фортрана есть статическая библиотека)
Праздный вопрос (праздный, потому как слинковали с фортраном): как в С делать "красивую" редукцию?
FORTRAN
SUBROUTINE CALCSTATS(A, N, &C++
max_real, min_real, &
max_imag, min_imag, &
max_abs, min_abs)
IMPLICIT NONE
integer(8), intent(in):: N
complex(4), intent(in):: A(N)
real(4), intent(inout):: max_real, min_real, &
max_imag, min_imag, &
min_abs, max_abs
real(4) my_abs
integer(8) i
!$OMP PARALLEL DO REDUCTION(MAX: max_real,max_imag,max_abs) REDUCTION(MIN: min_real,min_imag,min_abs) PRIVATE(my_abs)
do i = 1, N
max_real = MAX(REAL(A(i)), max_real)
min_real = MIN(REAL(A(i)), min_real)
max_imag = MAX(IMAG(A(i)), max_imag)
min_imag = MIN(IMAG(A(i)), min_imag)
my_abs = REAL(A(i)) * REAL(A(i)) + IMAG(A(i)) * IMAG(A(i))
max_abs = MAX(my_abs, max_abs)
min_abs = MIN(my_abs, min_abs)
end do
!$OMP END PARALLEL DO
max_abs = SQRT(max_abs)
min_abs = SQRT(min_abs)
end SUBROUTINE CALCSTATS
static void CalcStats(const std::complex<float> * A,
long long N,
float& max_real_val, float& min_real_val,
float& max_imag_val, float& min_imag_val,
float& min_abs_val, float& max_abs_val)
{
const int cores = omp_get_max_threads();
// Allocate
// It might crash with large # of cores :)
float * my_min_real_val = (float *)alloca(cores * sizeof(float));
float * my_max_real_val = (float *)alloca(cores * sizeof(float));
float * my_min_imag_val = (float *)alloca(cores * sizeof(float));
float * my_max_imag_val = (float *)alloca(cores * sizeof(float));
float * my_min_abs_val = (float *)alloca(cores * sizeof(float));
float * my_max_abs_val = (float *)alloca(cores * sizeof(float));
// Init
for (int core = 0; core < cores; core++)
{
my_min_real_val[core] = min_real_val;
my_max_real_val[core] = max_real_val;
my_min_imag_val[core] = min_imag_val;
my_max_imag_val[core] = max_imag_val;
my_min_abs_val[core] = min_abs_val;
my_max_abs_val[core] = max_abs_val;
}
#pragma omp parallel shared(cores, A)
{
// Init each core's TLS, to prevent cache line corruption
float my_min_real = min_real_val;
float my_max_real = max_real_val;
float my_min_imag = min_imag_val;
float my_max_imag = max_imag_val;
float my_min_abs = min_abs_val;
float my_max_abs = max_abs_val;
const int core = omp_get_thread_num();
long long start = (N / cores) * core;
long long end = (N / cores) * (core+1);
// Calculate on each core
for (long long i = start; i < end; i++)
{
float my_real = A[i].real();
float my_imag = A[i].imag();
float my_abs = (my_real * my_real + my_imag * my_imag);
my_min_real= min(my_real, my_min_real);
my_max_real= max(my_real, my_max_real);
my_min_imag= min(my_imag, my_min_imag);
my_max_imag= max(my_imag, my_max_imag);
my_min_abs = min(my_abs, my_min_abs);
my_max_abs = max(my_abs, my_max_abs);
}
// Reduction to common array
my_min_real_val[core] = my_min_real;
my_max_real_val[core] = my_max_real;
my_min_imag_val[core] = my_min_imag;
my_max_imag_val[core] = my_max_imag;
my_min_abs_val[core] = my_min_abs;
my_max_abs_val[core] = my_max_abs;
}
// Reduction to finals
for (int core = 0; core < cores; core++)
{
min_real_val = min(min_real_val, my_min_real_val[core]);
max_real_val = max(max_real_val, my_max_real_val[core]);
min_imag_val = min(min_imag_val, my_min_imag_val[core]);
max_imag_val = max(max_imag_val, my_max_imag_val[core]);
min_abs_val = min(min_abs_val, my_min_abs_val[core]);
max_abs_val = max(max_abs_val, my_max_abs_val[core]);
}
max_abs_val = sqrt(max_abs_val);
min_abs_val = sqrt(min_abs_val);
}
).
:
)
.
.
).
Всё-таки безобразие - жизненно важная конструкция не может быть реализована натуральным образом ! Только лишь потому, что для неё не придумали в своё время знака операции и приходится вызывать как функция. Я ведь правильно понимаю проблему ?
.
вставить выделенную цитату в окно ответа