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<vector<double>> (using std) 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

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 3

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 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

Answer 5

I don't recommend it, but not because of performance issues. It will be a little less performant than a traditional matrix, which are usually allocated as a big chunk of contiguous data that is indexed using a single pointer dereference and integer arithmetic . The reason for the performance hit is mostly caching differences, but once your matrix size gets big enough this effect will be amortized and if you use a special allocator for the inner vectors so that they are aligned to cache boundaries then this further mitigates the caching issue.

That by itself isn't reason enough not to do it, in my opinion. The reason for me is that it creates a lot of coding headaches. Here's a list of headaches this will cause long-term

Use of HPC libraries

If you want to use most HPC libraries you'll need to iterate over your vector and place all their data in a contiguous buffer, because most HPC libraries expect this explicit format. BLAS and LAPACK come to mind, but also the ubiquitous HPC library MPI would be much harder to use.

More potential for coding error

std::vector doesn't know anything about its entries. If you fill a std::vector with more std::vectors then it's entirely your job to make sure that they all have the same size, because remember that we want a matrix and matrices don't have variable number of rows (or columns). Thus you'll have to call all the correct constructors for every entry of your outer vector, and anybody else who uses your code must resist the temptation to use std::vector<T>::push_back() on any of the inner vectors, which would cause all following code to break. Of course you can disallow this if you write your class correctly, but it's much easier to enforce this simply with a big contiguous allocation.

HPC culture and expectations

HPC programmers simply expect low level data. If you give them a matrix there is an expectation that if they grabbed the pointer to the first element of the matrix and a pointer to the last element of the matrix, then all pointers in between these two are valid and point to elements of that same matrix. This is similar to my first point, but different because it may not be related so much to libraries but rather team members or anyone you share your code with.

Easier to reason about performance of lower level data

Dropping to the lowest level representation of your desired data structure makes your life easier in the long run for HPC. Using tools like perf and vtune will give you very low level performance counter measurements which you will try to combine with traditional profiling results in order to improve the performance of your code. If your data structure uses a lot of fancy containers it will be hard to understand that cache misses are coming from a problem with the container or an inefficiency in the algorithm itself. For more complicated code containers are necessary, but for matrix algebra they really aren't - you can get by with just 1 std::vector to store the data rather than n std::vectors, so go with that.

Answer 6

I also write a benchmark. For matrix of small size (<100*100), the performance is similar for vector< vector< double>> and wrapped 1D vector. For matrix of large size (~1000*1000), wrapped 1D vector is better. The Eigen matrix behaves worse. It is surprise to me that the Eigen is the worst.

#include <iostream>
#include <iomanip>
#include <fstream>
#include <sstream>
#include <algorithm>
#include <map>
#include <vector>
#include <string>
#include <cmath>
#include <numeric>
#include "time.h"
#include <chrono>
#include <cstdlib>
#include <Eigen/Dense>

using namespace std;
using namespace std::chrono;    // namespace for recording running time
using namespace Eigen;

int main()
{
    const int row = 1000;
    const int col = row;
    const int N = 1e8;

    // 2D vector
    auto start = high_resolution_clock::now();
    vector<vector<double>> vec_2D(row,vector<double>(col,0.));
    for (int i = 0; i < N; i++)
    {
        for (int i=0; i<row; i++)
        {
            for (int j=0; j<col; j++)
            {
                vec_2D[i][j] *= vec_2D[i][j];
            }
        }
    }
    auto stop = high_resolution_clock::now();
    auto duration = duration_cast<microseconds>(stop - start);
    cout << "2D vector: " << duration.count()/1e6 << " s" << endl;

    // 2D array
    start = high_resolution_clock::now();
    double array_2D[row][col];
    for (int i = 0; i < N; i++)
    {
        for (int i=0; i<row; i++)
        {
            for (int j=0; j<col; j++)
            {
                array_2D[i][j] *= array_2D[i][j];
            }
        }
    }
    stop = high_resolution_clock::now();
    duration = duration_cast<microseconds>(stop - start);
    cout << "2D array: " << duration.count() / 1e6 << " s" << endl;

    // wrapped 1D vector
    start = high_resolution_clock::now();
    vector<double> vec_1D(row*col, 0.);
    for (int i = 0; i < N; i++)
    {
        for (int i=0; i<row; i++)
        {
            for (int j=0; j<col; j++)
            {
                vec_1D[i*col+j] *= vec_1D[i*col+j];
            }
        }
    }
    stop = high_resolution_clock::now();
    duration = duration_cast<microseconds>(stop - start);
    cout << "1D vector: " << duration.count() / 1e6 << " s" << endl;

    // eigen 2D matrix
    start = high_resolution_clock::now();
    MatrixXd mat(row, col);
    for (int i = 0; i < N; i++)
    {
        for (int j=0; j<col; j++)
        {
            for (int i=0; i<row; i++)
            {
                mat(i,j) *= mat(i,j);
            }
        }
    }
    stop = high_resolution_clock::now();
    duration = duration_cast<microseconds>(stop - start);
    cout << "2D eigen matrix: " << duration.count() / 1e6 << " s" << endl;
}

Answer 7

As others have pointed out, don't try to do math with it or do anything performant.

That said, I've used this structure as a temporary when the code needs to assemble a 2-D array whose dimensions will be determined at runtime and after you've started storing data. For example, collecting vector outputs from some expensive process where it's not simple to compute exactly how many vectors you'll need to store at startup.

You could just concatenate all of your vector inputs into one buffer as they come in, but the code will be more durable and more readable if you use a vector<vector<T>>.