使用 vector<vector<double>> 形成高性能科学计算代码的矩阵类是个好主意吗?

计算科学 高性能计算 C++
2021-12-12 19:34:44

vector<vector<double>>使用(使用std)形成高性能科学计算代码的矩阵类是个好主意吗?

如果答案是否定的。为什么?谢谢

4个回答

这是一个坏主意,因为向量需要在空间中分配与矩阵中的行一样多的对象。分配是昂贵的,但主要是一个坏主意,因为矩阵的数据现在存在于分散在内存中的许多数组中,而不是全部在一个处理器缓存可以轻松访问它的地方。

这也是一种浪费的存储格式:std::vector 存储两个指针,一个指向数组的开头,一个指向结尾,因为数组的长度是灵活的。另一方面,要成为一个合适的矩阵,所有行的长度必须相同,因此只存储一次列数就足够了,而不是让每一行独立存储其长度。

除了 Wolfgang 提到的原因之外,如果您使用 a vector<vector<double> >,则每次要检索元素时都必须取消引用它两次,这比单个取消引用操作的计算成本更高。一种典型的方法是分配单个数组( avector<double>或 a double *)。我还看到人们通过围绕这个单个数组进行一些更直观的索引操作来向矩阵类添加语法糖,以减少调用正确索引所需的“心理开销”量。

不,使用免费的可用线性代数库之一。可以在这里找到关于不同库的讨论:推荐一个可用的、快速的 C++ 矩阵库?

真的有这么糟糕吗?

@Wolfgang:根据密集矩阵的大小,每行两个额外的指针可能可以忽略不计。关于分散数据,可以考虑使用自定义分配器,以确保向量位于连续内存中。只要内存没有被回收,即使是标准分配器,我们也会使用具有两个指针大小差距的连续内存。

@Geoff:如果您进行随机访问并且只使用一个数组,您仍然需要计算索引。可能不会更快。

所以让我们做一个小测试:

矢量矩阵.cc:

#include<vector>
#include<iostream>
#include<random>
#include <functional>
#include <sys/time.h>

int main()
{
  int N=1000;
  struct timeval start, end;

  std::cout<< "Checking differenz between last entry of previous row and first entry of this row"<<std::endl;
  std::vector<std::vector<double> > matrix(N, std::vector<double>(N, 0.0));
  for(std::size_t i=1; i<N;i++)
    std::cout<< "index "<<i<<": "<<&(matrix[i][0])-&(matrix[i-1][N-1])<<std::endl;
  std::cout<<&(matrix[0][N-1])<<" "<<&(matrix[1][0])<<std::endl;
  gettimeofday(&start, NULL);
  int k=0;

  for(int j=0; j<100; j++)
    for(std::size_t i=0; i<N;i++)
      for(std::size_t j=0; j<N;j++, k++)
        matrix[i][j]=matrix[i][j]*matrix[i][j];
  gettimeofday(&end, NULL);
  double seconds  = end.tv_sec  - start.tv_sec;
  double useconds = end.tv_usec - start.tv_usec;

  double mtime = ((seconds) * 1000 + useconds/1000.0) + 0.5;

  std::cout<<"calc took: "<<mtime<<" k="<<k<<std::endl;

  std::normal_distribution<double> normal_dist(0, 100);
  std::mt19937 engine; // Mersenne twister MT19937
  auto generator = std::bind(normal_dist, engine);
  for(std::size_t i=1; i<N;i++)
    for(std::size_t j=1; j<N;j++)
      matrix[i][j]=generator();
}

现在使用一个数组:

数组矩阵.cc

    #include<vector>
#include<iostream>
#include<random>
#include <functional>
#include <sys/time.h>

int main()
{
  int N=1000;
  struct timeval start, end;

  std::cout<< "Checking difference between last entry of previous row and first entry of this row"<<std::endl;
  double* matrix=new double[N*N];
  for(std::size_t i=1; i<N;i++)
    std::cout<< "index "<<i<<": "<<(matrix+(i*N))-(matrix+(i*N-1))<<std::endl;
  std::cout<<(matrix+N-1)<<" "<<(matrix+N)<<std::endl;

  int NN=N*N;
  int k=0;

  gettimeofday(&start, NULL);
  for(int j=0; j<100; j++)
    for(double* entry =matrix, *endEntry=entry+NN;
        entry!=endEntry;++entry, k++)
      *entry=(*entry)*(*entry);
  gettimeofday(&end, NULL);
  double seconds  = end.tv_sec  - start.tv_sec;
  double useconds = end.tv_usec - start.tv_usec;

  double mtime = ((seconds) * 1000 + useconds/1000.0) + 0.5;

  std::cout<<"calc took: "<<mtime<<" k="<<k<<std::endl;

  std::normal_distribution<double> normal_dist(0, 100);
  std::mt19937 engine; // Mersenne twister MT19937
  auto generator = std::bind(normal_dist, engine);
  for(std::size_t i=1; i<N*N;i++)
      matrix[i]=generator();
}

在我的系统上,现在有明显的赢家(编译器 gcc 4.7 和 -O3)

时间向量矩阵打印:

index 997: 3
index 998: 3
index 999: 3
0xc7fc68 0xc7fc80
calc took: 185.507 k=100000000

real    0m0.257s
user    0m0.244s
sys     0m0.008s

我们还看到,只要标准分配器不回收释放的内存,数据就是连续的。(当然,在一些解除分配之后,不能保证这一点。)

时间数组矩阵打印:

index 997: 1
index 998: 1
index 999: 1
0x7ff41f208f48 0x7ff41f208f50
calc took: 187.349 k=100000000

real    0m0.257s
user    0m0.248s
sys     0m0.004s