【发布时间】：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 回复：stackoverflow.com/questions/30838959/…。请注意，我有点“不建议”使用内在函数（因为可移植性）进行 simdization + 这篇文章确实不是关于您的问题陈述，但这仍然是考虑 2 级并行性以及如何处理的好方法记忆墙。
标签： c++ matrix graph parallel-processing