通过内在函数将 Doubles 累积到 Bins 中

Question

我有一个观测值向量和一个等长的偏移量向量，将观测值分配给一组 bin。每个 bin 的值应该是分配给该 bin 的所有观察值的总和，我想知道是否有矢量化方法来进行归约。

下面是一个简单的实现：

const int N_OBS = 100`000`000;
const int N_BINS = 16;
double obs[N_OBS];    // Observations
int8_t offsets[N_OBS];
double acc[N_BINS] = {0};

for (int i = 0; i < N_OBS; ++i) {
  acc[offsets[i]] += obs[i]; // accumulate obs value into its assigned bin
}

是否可以使用 simd/avx 内在函数？类似于上面的内容将被运行数百万次。我查看了 scatter/gather 方法，但似乎无法找到完成它的好方法。

Answer 1

现代 CPU 非常好运行你的幼稚版本。在 AMD Zen3 上，输入 100M 随机数需要 48 毫秒，即 18 GB/sec RAM 读取带宽。这大约是我计算机 (dual-channel DDR4-3200) 硬带宽限制的 35%。

恐怕没有 SIMD 可以提供帮助。尽管如此，我得到的最好的版本如下。使用 OpenMP 支持进行编译，开关取决于您的 C++ 编译器。

void computeHistogramScalarOmp( const double* rsi, const int8_t* indices, size_t length, double* rdi )
{
    // Count of OpenMP threads = CPU cores to use
    constexpr int ompThreadsCount = 4;

    // Use independent set of accumulators per thread, otherwise concurrency gonna corrupt data.
    // Aligning by 64 = cache line, we want to assign cache lines to CPU cores, sharing them is extremely expensive
    alignas( 64 ) double accumulators[ 16 * ompThreadsCount ];
    memset( &accumulators, 0, sizeof( accumulators ) );

    // Minimize OMP overhead by dispatching very few large tasks
#pragma omp parallel for schedule(static, 1)
    for( int i = 0; i < ompThreadsCount; i++ )
    {
        // Grab a slice of the output buffer
        double* const acc = &accumulators[ i * 16 ];

        // Compute a slice of the source data for this thread
        const size_t first = i * length / ompThreadsCount;
        const size_t last = ( i + 1 ) * length / ompThreadsCount;

        // Accumulate into thread-local portion of the buffer
        for( size_t i = first; i < last; i++ )
        {
            const int8_t idx = indices[ i ];
            acc[ idx ] += rsi[ i ];
        }
    }

    // Reduce 16*N scalars to 16 with a few AVX instructions
    for( int i = 0; i < 16; i += 4 )
    {
        __m256d v = _mm256_load_pd( &accumulators[ i ] );
        for( int j = 1; j < ompThreadsCount; j++ )
        {
            __m256d v2 = _mm256_load_pd( &accumulators[ i + j * 16 ] );
            v = _mm256_add_pd( v, v2 );
        }
        _mm256_storeu_pd( rdi + i, v );
    }
}

以上版本导致 20.5 毫秒的时间，转化为 88% 的 RAM 带宽限制。

P.S。我不知道为什么这里的最佳线程数是 4，我在 CPU 中有 8 cores/16 个线程。较低和较高的值都会降低带宽。常量大概是CPU-specific.

Answer 2

如果 offsets 确实没有改变数千次（可能甚至几十次），那么“转置”它们可能是值得的，即存储所有需要添加到 [= 的索引12=]，然后是需要添加到 acc[1] 的所有索引，等等

本质上，您最初所做的是 sparse-matrix 乘以 dense-vector 的乘积，矩阵采用 compressed-column-storage 格式（没有显式存储 1 值）。

如所示，如果矩阵存储在compressed-row-storage中，稀疏GEMV产品通常会更快（即使没有AVX2的收集指令，您也不需要每次都加载和存储累加值时间）。

未经测试的示例实现：

using sparse_matrix = std::vector<std::vector<int> >;

// call this once:
sparse_matrix transpose(uint8_t const* offsets, int n_bins, int n_obs){
    sparse_matrix res;
    res.resize(n_bins);

    // count entries for each bin: 
    for(int i=0; i<n_obs; ++i) {
        // assert(offsets[i] < n_bins);
        res[offsets[i]].push_back(i);
    }

    return res;
}

void accumulate(double acc[], sparse_matrix const& indexes, double const* obs){
    for(std::size_t row=0; row<indexes.size(); ++row) {
        double sum = 0;
        for(int col : indexes[row]) {
            // you can manually vectorize this using _mm256_i32gather_pd,
            // but clang/gcc should autovectorize this with -ffast-math -O3 -march=native
            sum += obs[col];
        }
        acc[row] = sum;
    }
}

通过内在函数将 Doubles 累积到 Bins 中

Accumulating Doubles Into Bins via intrinsics

c++

simd

avx

avx2