Предполагая, что вы используете современный процессор Intel/AMD (с AVX), и вы хотите сохранить тот же алгоритм, вы можете попробовать следующий код. Он использует AVX и OpenMP для распараллеливания. Скомпилируйте с помощью GCC foo.c -mavx -fopenmp -O3. Если вы не хотите использовать OpenMP, просто закомментируйте два оператора #pragma.
Скорость будет зависеть от размеров массива и размеров кеша. Для массивов, которые вписываются в кеш L1, вы можете ожидать около 6-кратного ускорения (вы должны отключить OpenMP, а затем из-за его накладных расходов). Повышение будет продолжаться с каждым уровнем кеша. Когда он добирается до системной памяти, он все равно получает импульс, хотя (работа над удвоением 10M (2 * 80 МБ) по-прежнему на 70% выше по сравнению с моей двухжильной мостовой системой плюща).
#include <stdio.h>
#include <stdlib.h>
#include <immintrin.h>
#include <omp.h>
double foo_avx_omp(const double *descr1, const double *descr2, const int d) {
double diff, dsq = 0;
int i;
int range;
__m256d diff4_v1, diff4_v2, dsq4_v1, dsq4_v2, t1, l1, l2, l3, l4;
__m128d t2, t3;
range = (d-1) & -8; //round down to multiple of 8
#pragma omp parallel private(i,l1,l2,l3,l4,t1,t2,t3,dsq4_v1,dsq4_v2,diff4_v1,diff4_v2) \
reduction(+:dsq)
{
dsq4_v1 = _mm256_set1_pd(0.0);
dsq4_v2 = _mm256_set1_pd(0.0); //two sums to unroll the loop once
#pragma omp for
for(i=0; i<(range/8); i++) {
//load one cache line of descr1
l1 = _mm256_load_pd(&descr1[8*i]);
l3 = _mm256_load_pd(&descr1[8*i+4]);
//load one cache line of descr2
l2 = _mm256_load_pd(&descr2[8*i]);
l4 = _mm256_load_pd(&descr2[8*i+4]);
diff4_v1 = _mm256_sub_pd(l1, l2);
diff4_v2 = _mm256_sub_pd(l3, l4);
dsq4_v1 = _mm256_add_pd(dsq4_v1, _mm256_mul_pd(diff4_v1, diff4_v1));
dsq4_v2 = _mm256_add_pd(dsq4_v2, _mm256_mul_pd(diff4_v2, diff4_v2));
}
dsq4_v1 = _mm256_add_pd(dsq4_v1, dsq4_v2);
t1 = _mm256_hadd_pd(dsq4_v1,dsq4_v1);
t2 = _mm256_extractf128_pd(t1,1);
t3 = _mm_add_sd(_mm256_castpd256_pd128(t1),t2);
dsq += _mm_cvtsd_f64(t3);
}
//finish remaining elements if d was not a multiple of 8
for (i=range; i < d; ++i) {
diff = descr1[i] - descr2[i];
dsq += diff * diff;
}
return dsq;
}
double foo(double *descr1, double *descr2, int d) {
double diff, dsq = 0;
int i;
for (i = 0; i < d; ++i)
{
diff = descr1[i] - descr2[i];
dsq += diff * diff;
}
return dsq;
}
int main(void)
{
double result1, result2, result3, dtime;
double *descr1, *descr2;
const int n = 2000000;
int i;
int repeat = 1000;
descr1 = _mm_malloc(sizeof(double)*n, 64); //align to a cache line
descr2 = _mm_malloc(sizeof(double)*n, 64); //align to a cache line
for(i=0; i<n; i++) {
descr1[i] = 1.0*rand()/RAND_MAX;
descr2[i] = 1.0*rand()/RAND_MAX;
}
dtime = omp_get_wtime();
for(i=0; i<repeat; i++) {
result1 = foo(descr1, descr2, n);
}
dtime = omp_get_wtime() - dtime;
printf("foo %f, time %f\n", result1, dtime);
dtime = omp_get_wtime();
for(i=0; i<repeat; i++) {
result1 = foo_avx_omp(descr1, descr2, n);
}
dtime = omp_get_wtime() - dtime;
printf("foo_avx_omp %f, time %f\n", result1, dtime);
return 0;
}
Можете ли вы переключиться на целые числа? Какой диапазон значений будет охватываться содержимым 'descr1/2'? – alk
Вы не инициализировали d! –
Если зависимости позволяют запускать ядро на графическом процессоре – jev