C语言 改进递归哈达马尔变换

Question

我有以下代码来计算 Hadamard 变换。现在，hadamard函数是我的程序的瓶颈。你认为有任何加速的潜力吗？也许使用 AVX2 指令？典型的输入大小约为 512 或 1024。

贝斯特，汤姆

·

#include <stdio.h>
void hadamard(double *p, size_t len) {
    double tmp = 0.0;
    if(len == 2) {
        tmp = p[0];
        p[0] = tmp + p[1];
        p[1] = tmp - p[1];
    } else {
        hadamard(p, len/2);
        hadamard(p+len/2, len/2);
        for(int i = 0; i < len/2; i++) {
           tmp = p[i];
           p[i] = tmp + p[i+len/2];
           p[i+len/2] = tmp - p[i+len/2];
       }
   }
}
int main(int argc, char* argv[]) {
        double a[] = {1.0, 2.0, 3.0, 4.0};
        hadamard(a, 4);
}

Answer 1

这是一个概念验证实现，基于快速沃尔什-哈达玛变换，具有优化的首遍。 使用 clang-3.3 或更高版本编译良好，使用 clang-4.0 或更高版本提供不错的结果（需要 -O2 才能正确内联相关函数）。如果你没有FMA，你需要在hadamard4中用-0.0 hada2_的下两个元素进行异化（并使用普通_mm256_add_pd）。

我检查的所有 gcc 版本都需要通过手动加载/存储内部函数替换memcpy以获得类似的结果。

此外，我把处理案件len<16作为练习。并且可能值得实施一个hadamard32，也许是一个类似于hadamard16 hadamard64，以更好地使用可用的寄存器并减少内存访问。C++这可以通过递归模板实现来完成。

#include <immintrin.h> // avx+fma
#include <assert.h> // assert
#include <string.h> // memcpy
inline __m256d hadamard4(__m256d x0123)
{
    __m256d x1032 = _mm256_permute_pd(x0123, 5);             // [x1, x0, x3, x2]
    __m256d hada2 = _mm256_addsub_pd(x1032,x0123);           // [x0+x1, x0-x1, x2+x3, x2-x3]
    __m256d hada2_= _mm256_permute2f128_pd(hada2, hada2, 1); // [x2+x3, x2-x3, x0+x1, x0-x1]
    // if no FMA is available, this can be done with xoring and adding:
    __m256d res   = _mm256_fmadd_pd(hada2_, _mm256_set_pd(1.0, 1.0, -1.0, -1.0), hada2);
    return res;
}
inline void hadamard8(__m256d data[2])
{
    __m256d a = hadamard4(data[0]);
    __m256d b = hadamard4(data[1]);
    data[0] = _mm256_add_pd(a,b);
    data[1] = _mm256_sub_pd(a,b);
}
inline void hadamard16(__m256d data[4])
{
    hadamard8(data+0);
    hadamard8(data+2);
    for(int i=0; i<2; ++i)
    {
        __m256d tmp = data[i];
        data[i]  = _mm256_add_pd(tmp, data[i+2]);
        data[i+2]= _mm256_sub_pd(tmp, data[i+2]);
    }
}
void hadamard(double* p, size_t len)
{
    assert((len&(len-1))==0); // len must be power of 2
    assert(len>=16); // TODO implement fallback for smaller sizes ...
    // first pass: hadamard of 16 values each
    for(size_t i=0; i<len; i+=16)
    {
        __m256d data[4];
        memcpy(data, p+i, sizeof(data)); // should get optimized to 4x vmovupd
        hadamard16(data);
        memcpy(p+i, data, sizeof(data)); // should get optimized to 4x vmovupd
    }
    for(size_t h=32; h<len; h*=2)
    {
        for(size_t i=0; i<len; i+=2*h)
        {
            for(size_t j=i; j<i+h; j+=4)
            {
                __m256d x = _mm256_loadu_pd(p+j);
                __m256d y = _mm256_loadu_pd(p+j+h);
                _mm256_storeu_pd(p+j,   _mm256_add_pd(x,y));
                _mm256_storeu_pd(p+j+h, _mm256_sub_pd(x,y));
            }
        }
    }
}

基准测试也是一种练习;-）

免责声明：我没有测试这个。看着它，我可能在_mm256_fmadd_pd指令中混淆了hada2_和hada2......