double diff, dsq = 0;
double *descr1, *descr2;
int i, d;
for (i = 0; i < d; ++i)
{
diff = descr1[i] - descr2[i];
dsq += diff * diff;
}
return dsq;
我想优化这部分代码，这在我的程序中花费了大部分时间。如果这种双重乘法以优化的方式执行，我的程序可以运行得非常快。是否有其他乘法方法而不是使用 * 运算符导致程序运行得更快？多谢。
user2758590
问问题
363 次
9 回答
3
这绝对是Duff's Device 的一个案例。
这是我的实现，基于 Duff 的设备。
（注意：仅经过轻微测试......必须在调试器中逐步完成以确保正确的行为）
void fnc(void)
{
double dsq = 0.0;
double diff[8] = {0.0};
double descr1[115];
double descr2[115];
double* pD1 = descr1;
double* pD2 = descr2;
int d = 115;
//Fill with random data for testing
for(int i=0; i<d; ++i)
{
descr1[i] = (double)rand() / (double)rand();
descr2[i] = (double)rand() / (double)rand();
}
// Duff's Device: Step through this in a debugger, its AMAZING.
int c = (d + 7) / 8;
switch(d % 8) {
case 0: do { diff[0] = *pD1++ - *pD2++; diff[0] *= diff[0];
case 7: diff[7] = *pD1++ - *pD2++; diff[7] *= diff[7];
case 6: diff[6] = *pD1++ - *pD2++; diff[6] *= diff[6];
case 5: diff[5] = *pD1++ - *pD2++; diff[5] *= diff[5];
case 4: diff[4] = *pD1++ - *pD2++; diff[4] *= diff[4];
case 3: diff[3] = *pD1++ - *pD2++; diff[3] *= diff[3];
case 2: diff[2] = *pD1++ - *pD2++; diff[2] *= diff[2];
case 1: diff[1] = *pD1++ - *pD2++; diff[1] *= diff[1];
dsq += diff[0] + diff[1] + diff[2] + diff[3] + diff[4] + diff[5] + diff[6] + diff[7];
} while(--c > 0);
}
}
解释
正如其他人所说，您几乎无法优化浮点运算。但是，在您的原始代码中，程序花费了大量时间检查
i.
执行步骤大致为：
Is i < d? ==> Yes
Do some math.
Is i < d? ==> Yes
Do some math.
Is i < d? ==> Yes
Do some math.
Is i < d? ==> Yes
Do some math.
您可以看到其他所有步骤都在检查i。
使用 Duff 的设备，您可以在检查计数器之前进行八次c操作（在本例中）。
现在执行步骤大致是：
Is c > 0? ==> Yes
Do some math.
Do some math.
Do some math.
Do some math.
Do some math.
Do some math.
Do some math.
Do some math.
Is c > 0? ==> Yes
Do some math.
Do some math.
Do some math.
Do some math.
Do some math.
Do some math.
Do some math.
Do some math.
Is c > 0? ==> Yes
[...]
换句话说，您实际完成工作所花费的 CPU 大约是 8 倍，而检查计数器值的时间要少得多。这是一个巨大的胜利。
我怀疑您甚至可以将循环进一步展开到 16 或 32 次操作，以获得更大的胜利。这实际上取决于d代码中的可能值。
请测试并分析此代码，并让我知道它对您的效果如何。
我有一种强烈的感觉，这将是一个很大的进步。
于 2013-09-08T12:45:02.953 回答
2
您可以通过其严格的别名规则帮助编译器：
double calc_ssq(double *restrict descr1, double *restrict descr2, size_t count)
{
double ssq;
ssq = 0.0;
for ( ;count; count--) {
double diff;
diff = *descr1++ - *descr2++;
ssq += diff * diff;
}
return ssq;
}
于 2013-09-08T11:49:06.353 回答
1
假设您使用的是现代 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;
}
于 2013-09-08T17:37:50.467 回答
1
总而言之，您有一个非常紧密的循环来访问大量数据。循环展开可能有助于隐藏延迟，但在现代硬件上，像这样的循环受限于内存带宽，而不是计算能力。
因此，您对优化的唯一真正希望是：a) 使用数组float而不是数组double来将从内存加载的数据量减少一半，以及 b) 尽可能避免调用此代码。
以下是一些数字：
您的内部循环中有三个双算术指令，大约是 6 个周期。这些需要 16 个字节的数据。在 3 GHz 处理器上，内存带宽为 8 GB/s。DDR3-1066 模块提供 8.5 GB/s。因此，即使您使用 SSE 之类的东西，也不会变得更快，除非您切换到使用float.
于 2013-09-08T13:34:47.200 回答
1
看起来您正在计算两个向量的均方误差。
使用BLAS，您将能够利用手动优化的代码，其效率远超我们任何人编写的。
于 2013-09-09T16:09:32.807 回答
1
如果您的计算并不真的需要双精度，您可以尝试将它们转换为单精度，然后相乘。
我想，在 32 位处理器的情况下，单精度乘法将比双精度乘法更快，因为常规float只需要一个处理器寄存器并且double需要两个。
我不确定铸造不会“吃掉”所有速度改进，你会从单精度乘法中获得。
于 2013-09-08T11:49:45.430 回答
1
如果 d 可以被 2 整除，我会尝试这样的事情：
for(i=0;i<d;i+=2)
{
diff0 = descr1[i] - descr2[i];
diff1 = descr1[i+1] - descr2[i+1];
dsq += diff0 * diff0 + diff1 * diff1;
}
这将暗示优化器可以交错六个操作。即使 d 是奇数，您也可以在每个向量的末尾附加一个 0.0 值（给出偶数个值），因为考虑到所涉及的操作，这对结果没有影响。
下一步可能是将向量附加到可被 4 整除，在 i+=4 迭代之前进行四次减法、四次乘法和四次加法；
能被 8 整除的向量可以完全适合 64 的高速缓存行大小。
浮点乘法只需要一个或两个时钟周期即可完成，加法和减法也是如此（根据 Agner Fog 的说法）。因此，对于您的示例，减少迭代开销应该可以加快速度。
于 2013-09-08T12:23:40.610 回答
0
警告：以下未经测试的代码。
如果您的硬件和编译器都支持它们，您可能希望使用向量并行化某些操作。我过去在 GCC 4.6.x 编译器（x86-64 Ubuntu 机器）上使用过类似于以下内容的东西。如果使用不同的编译器/体系结构，某些语法可能会略有不同或可能会有所不同。但是，我希望它足够接近以使您更接近您的目标。
typedef double v2d_t __attribute__((vector_size (16)));
typedef union {
v2d_t vector;
double d[2];
} v2d_u;
v2d_u vdsq = (v2d_t) {0.0, 0.0}; /* sum of square of differences */
v2d_u vdiff; /* difference */
v2d_t * vdescr1; /* pointer to array of aligned vector of doubles */
v2d_t * vdescr2; /* pointer to array of aligned vector of doubles */
int i; /* index into array of aligned vector of doubles */
int d; /* # of elements in array */
/*
* ...
* Assuming that <d> is getting initialized appropriately somewhere
* ...
*/
for (i = 0; i < d; i++) {
vdiff.vector = vdescr1[i] - vdescr2[i];
vdsq.vector += vdiff.vector * vdiff.vector;
}
return vdsq.d[0] + vdsq.d[1];
上面可能可以进行更多调整以获得更好的性能。也许一些循环展开。或者，如果您可以利用 256 位向量（例如某些 x86 处理器上的 YMMx）而不是此示例使用的 128 位向量，那也可能会加快速度（需要对代码进行一些调整）。
希望这可以帮助。
于 2013-09-08T13:50:43.257 回答
0
descr1将和的两个替身descr2放入一个结构中，使它们在内存中彼此相邻。这将使缓存的使用和对内存的访问更好。
也register用于diff和dsq
于 2013-09-08T10:56:49.120 回答