Is it a good idea to use vector<vector<double>> to form a matrix class for high performance scientific computing code?

Question

Is it a good idea to use vector> to form a matrix class for high performance scientific computing code?

If the answer is no. Why? Thanks

Answer 1

It's a bad idea because vector needs to allocate as many objects in space as there are rows in your matrix. Allocation is expensive, but primarily it is a bad idea because the data of your matrix now exists in a number of arrays scattered around memory, rather than all in one place where the processor cache can easily access it.

It's also a wasteful storage format: std::vector stores two pointers, one to the beginning of the array and one to the end because the length of the array is flexible. On the other hand, for this to be a proper matrix, the lengths of all rows must be the same and so it would be sufficient to store the number of columns only once, rather than letting each row store its length independently.

Answer 2

No, use one of the free available linear algebra libraries. A discussion about different libraries can be found here: Recommendations for a usable, fast C++ matrix library?

Answer 3

In addition to the reasons Wolfgang mentioned, if you use a vector<vector<double> >, you'll have to dereference it twice every time you want to retrieve an element, which is more computationally costly than a single dereferencing operation. One typical approach is to allocate a single array (a vector<double> or a double *) instead. I've also seen people add syntactic sugar to matrix classes by wrapping around this single array some more intuitive indexing operations, to reduce the amount of "mental overhead" needed to invoke the proper indices.

Answer 4

Is it really such a bad thing?

@Wolfgang: Depending on the size of the dense matrix, two additional pointer per row might be negligable. Concerning scattered data one could think of using a custom allocator that makes sure that the vectors are in contiguous memory. As long as memory is not recycled even the standard allocator will us contiguous memory with a two pointer size gap.

@Geoff: If you are doing random access and use just one array you still have to calculate the index. Might not be faster.

So let us do a small test:

vectormatrix.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();
}

And now using one array:

arraymatrix.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();
}

On my system there is now clear winner (Compiler gcc 4.7 with -O3)

time vectormatrix prints:

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

We also see, that as long as the standard allocator does not recycle freed memory, the data is contiguous. (Of course after some deallocations there is no guarantee for this.)

time arraymatrix prints:

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