如何加快稀疏数组加法

Posted 2023-02-16

【中文标题】如何加快稀疏数组加法

【英文标题】：How to speed up sparse array addition

【发布时间】：2015-06-11 17:44:41

【问题描述】：

出于数据模拟的目的，我正在寻找一种有效的方法来对稀疏矩阵进行加权求和。 基本上我有一个 Nx x Ny x Nz 双值的数据立方体，其中 Nx 和 Ny 大约是 4000，而 Nz 是几百万。所有 Nx x Ny 子矩阵都非常稀疏（大约 40 个整数的数据块）。 我现在想通过将所有矩阵相加并加权它们来减少 Z 方向的数据立方体。该过程如图所示。对于我的模拟，所有矩阵都保持不变，只有权重会改变并生成不同的 Nx x Ny 数据集。

这是我尝试过的：C++ 中稀疏矩阵的简单实现，以及简单的求和。

#ifndef SPARSEARRAY3D_H
#define SPARSEARRAY3D_H

#include <vector>

struct data
    unsigned short int x;
    unsigned short int y;
    int z;
    double value;
;

class sparsearray3d

public:
    sparsearray3d();
    void createRandomData(int Nx, int Ny, int Nz, int bwidthX, int bwidthY);
    void sumData();
    int Nx,Ny,Nz;

    std::vector<data> dd;
    std::vector<std::vector<double> > image;
    std::vector<double> weights;
;

#endif // SPARSEARRAY3D_H

sparsearray3d.cpp

#include "sparsearray3d.h"
#include <stdlib.h>     /* srand, rand */
#include <stdio.h>      /* printf, scanf, puts, NULL */


sparsearray3d::sparsearray3d()

    this->Nx = 0;
    this->Ny = 0;
    this->Nz = 0;


void sparsearray3d::createRandomData(int Nx, int Ny, int Nz, int bwidthX = 5, int bwidthY = 5)

    // create random data
    this->weights.resize(Nz);

    this->image.resize( Nx , std::vector<double>( Ny , 0. ) );

    this->Nx = Nx;
    this->Ny = Ny;
    this->Nz = Nz;

    for(int i=0; i<Nz; ++i)
    
        int x0 = rand() % (Nx-bwidthX);
        int y0 = rand() % (Ny-bwidthY);

        this->weights.push_back((double) rand() / (RAND_MAX));

        for(int j=0; j<bwidthX; ++j)
        
            for(int k=0; k<bwidthY; ++k)
            
                this->dd.push_back(x0+j,y0+k,i,((double) rand() / (RAND_MAX)));
            
        
    
    printf("Vector size: %4.2f GB \n", this->dd.size()*sizeof(data) * 1E-9);



void sparsearray3d::sumData()

    std::vector<data>::iterator it;
    #pragma omp parallel for
    for(it = this->dd.begin(); it < this->dd.end(); ++it)
    
        this->image[it->y][it->x] += it->value * this->weights[it->z];
    

main.cpp

#include <iostream>
#include "sparsearray3d.h"
#include <sys/time.h>

using namespace std;

int main()


struct timeval start, end;

sparsearray3d sa;
gettimeofday(&start, NULL);
sa.createRandomData(4096, 4096, 2000000, 4, 16);
gettimeofday(&end, NULL);
double delta = ((end.tv_sec  - start.tv_sec) * 1000000u +
         end.tv_usec - start.tv_usec) / 1.e6;

cout << "random array generation: " << delta << endl;
gettimeofday(&start, NULL);
sa.sumData();
gettimeofday(&end, NULL);
delta = ((end.tv_sec  - start.tv_sec) * 1000000u +
         end.tv_usec - start.tv_usec) / 1.e6;
cout << "array addition: " << delta << endl;
return 0;

这已经做得很好了，上面的例子在这里运行大约 0.6 秒。 我想知道的第一件事是，为什么 #pragma omp parallel for 尽管使用了 4 个 CPU，但速度却只有 2 倍左右。

这个问题似乎非常适合大规模并行化。 Cuda / OpenCL 可以在这里提供帮助吗？但是，我在某处读到 Cuda/OpenCL 的矩阵加法不是很有效。 （虽然我没有可用的 NVIDIA 卡）。 或者，我阅读了一些关于图形及其与矩阵的关系的内容。这个问题可以用一些图算法来解决吗？

编辑： 我试图给 Eigen 一个机会；但是，我未能创建大量矩阵。下面的代码比我的代码需要更多的内存（N~20000000 失败，因为我的内存用完了）。不确定我做得对，但这就是我从 eigen 文档中理解的方式。

#include <vector>
#include <eigen3/Eigen/Sparse>

int main()

int N=100000;
std::vector<Eigen::SparseMatrix<double> > data;

data.resize(N);

for (int i=0; i<N; ++i)

    data[i].resize(4096,4096);
    data[i].reserve(4*16);


return 0;

另外，用以下方式对稀疏矩阵求和比我的代码慢得多：

Eigen::SparseMatrix<double> sum(4096,4096) ;
sum.reserve(4096*4096);
for(int i=0; i<N; ++i)
    sum+=data[i];

【问题讨论】：

你为什么不使用库（可能是 Boost、Eigen 等）？

查看编辑（本征耗尽了我的记忆，而且似乎很慢，至少在我使用它的（很可能是错误的）方式上）

也许您可以创建项目列表矩阵？这些矩阵的每个元素只是列表（或数组或其他东西），其中包含每个稀疏矩阵的 i,j 项。这些方法允许您对每个项目进行独立于其他项目的总和，并且 GPU 上的计算可以为您提供很好的加速。您还调整了 dd "this->dd.resize" 的大小，然后添加了新元素。确定吗？ PS对不起我的英语不好

谢谢，我更正了代码。我一开始没有使用 pushback，这就是为什么我仍然有 dd.resize 的原因。我会试试你的方法。

除了我的一般性和有意抽象的答案之外，您还可以在这里查看一个（在许多类似的）漂亮的 Z Boson 回复：***.com/questions/30838959/…。请注意，我有点“不建议”使用内在函数（因为可移植性）进行 simdization + 这篇文章确实不是关于您的问题陈述，但这仍然是考虑 2 级并行性以及如何处理的好方法记忆墙。

【参考方案1】：

您处理的是线性代数矩阵向量乘法的相对一般情况（要搜索的一些相关计算机科学缩写是“DGEMM”或“SpMV”）。因此，您的第一个选择是尝试高度优化的并行线性代数库，例如英特尔 MKL，另请参阅： Parallel linear algebra for multicore system

其次，如果您想自己优化和并行化您的算法，那么您可能首先需要学习一些相关文章（例如： - http://cscads.rice.edu/publications/pdfs/Williams-OptMultiCore-SC07.pdf - http://pcl.intel-research.net/publications/ics26-liuPS.pdf )

第三，如果您不想或没有时间浏览书籍、文章或图书馆，但想自己从头开始尝试一切，则需要考虑几件事：

您所拥有的实际上是细粒度 并行性（每次迭代都太快 == 太小），因此标准化的“每次迭代”线程调度开销可能相对难以“摊销”。 对于细粒度并行，最好尝试 SIMD 并行化（尤其是因为您有可能在此处使用“FMA==fused multiply add”）。在 OpenMP4.x（由 ICC 和 GCC4.9+ 支持）中，您有 #pragma omp simd ，它非常适合减少 simd。 但正确的选择实际上是简化/逆向工程您的示例，并使用 x 的 for 循环 (#pragma omp for) 和 y 的下一个循环 (#pragma omp simd) 使其显式循环。然后您将拥有 2 级并行度，这将使您的状态更好。 一旦您成功完成上述所有操作，您很快就会受到内存缓存/延迟、内存带宽或不规则访问模式的限制 - 类似于 内存墙 的 3 个面.在您当前的实现中，我希望您受到带宽的限制，但是真正优化的稀疏矩阵矩阵产品往往受到给定 3 个面的各种组合的限制。这是需要阅读一些文章的地方，考虑到“循环阻塞”、“数组结构”、“预取和流式传输”主题将变得非常相关。

最后（也许是第一项），我不认为您在这里真的有“稀疏矩阵”感知数据结构，因为您并没有真正尝试避免“零”（它们仍然位于您的矩阵中-向量）。我不知道你是否已经采用了你的数据结构（所以你的例子只是简化了真实的代码..）。

并行线性代数是计算机科学中的一个大课题；每年都有一些专家提出新的情报。每个下一代 CPU 和有时协处理器都经过更好的优化，以加速并行化线性代数内核。所以有很多东西可以探索。