假设您使用的是现代 Intel/AMD 处理器（带有 AVX）并且您希望保持相同的算法，您可以尝试以下代码。它使用 AVX 和 OpenMP 进行并行化。使用GCC foo.c -mavx -fopenmp -O3 编译。如果您不想使用 OpenMP，只需将两个 #pragma 语句注释掉即可。
速度将取决于数组大小和缓存大小。对于适合 L1 缓存的阵列，您可以预期大约 6 倍的加速（由于其开销，您应该禁用 OpenMP）。提升将随着每个缓存级别而不断下降。当它到达系统内存时，它仍然会得到提升（在我的两个核心常春藤桥系统上运行超过 10M 双倍 (2*80MB) 仍然快 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;
}