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

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

-
-1 Of course it's a bad idea. You won't be able to use blas, lapack or any other existing matrix library with such a storage format. In addition, you introduce inefficiencies by data non-localty and indirection –  Thomas Klimpel Aug 29 '12 at 6:58
@Thomas Does that really warrant a downvote? –  akid Aug 29 '12 at 8:12
Don't downvote. It's a legitimate question even if it's a misguided idea. –  Wolfgang Bangerth Aug 29 '12 at 11:02
std::vector is not a distributed vector so you won't be able to do much parallel computing with it (except for shared memory machines), use Petsc or Trilinos instead. Furthermore one usually deals with sparse matrices and you would be storing full dense Matrices. For playing with sparse matrices you could use a std::vector<std::map> but again, this would not perform very good, see @WolfgangBangerth post below. –  gnzlbg Aug 29 '12 at 12:13
thanks, I know it is not a good idea. Just want to know the exact reason. –  cfdgeek Aug 29 '12 at 20:50

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.

-

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?

-

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.

-

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

-
You write "On my system there is now clear winner" - did you mean no clear winner? –  akid Sep 1 '12 at 12:17
-1 Understanding the performance of hpc code can be nontrivial. In your case, the size of the matrix simply exceed the cache size, so that you are just measuring the memory bandwidth of your system. If I change N to 200 and increase the number of iterations to 1000, I get "calc took: 65" vs "calc took: 36". If I further replace a=a*a by a+=a1*a2 to make it more realistic, I get "calc took: 176" vs "calc took: 84". So it looks like you can loose a factor two in performance by using a vector of vectors instead of a matrix. Real life will be more complicated, but it's still a bad idea. –  Thomas Klimpel Sep 2 '12 at 8:53
yeah but try using std::vectors with MPI, C wins hands down –  pyCthon Sep 17 '12 at 4:31