2

I am considering implementing an array like container, and I'm not sure whether to use a gsl::gsl_vector or a std::vector. The container needs to be space efficient but also very fast at calling values. The container will be referenced constantly in the main program to, for example, input values into tensor functions, amongst other things. I am calling from the container literally billions of times.

Here are the pros and cons that I've considered so far: The gsl_vector is convenient because it will allow me to use the gsl BLAS libraries on occasion, and the gsl_vector_get(...) call is very efficient. On the other hand, I am able to get almost the same speed of calls using stl iterators, and the stl vector has an interface that I find quite natural.

Are there any memory overheads/efficiency issues I should be aware of that I've overlook in the above code?

Also, I am using a std::vector<gsl_vector*> implementation at the moment, and an iterator to traverse the std::vector. Would it be smarter to use a gsl_matrix here? The idea would be to use gsl_vector_views to get the right vector, rather than the iterator. Would this be more efficient?

Plamen
  • 650
  • 1
  • 8
  • 27

1 Answers1

2

On one hand, it is true that with gsl_vector you can use gsl BLAS which is a big advantage. On the other hand, it is also true that gsl interface is quite cumbersome for a c++ programmer. So, neither solutions are truly satisfactory. However, I strongly prefer the use of gsl_matrix because

(i) with some effort you can write a small wrapper class that ameliorate the cumbersome C interface of gsl_matrix (it is much harder to deal with the lack of BLAS library in std::vector).

(ii) gsl_matrix is just a wrapper to an one dimensional continuous array where m(i,j) = array[i*N + j] for square matrix (even if the matrix is not square gsl_matrix still implements it as one dimensional array). In std::vector<gsl_vector*>, you will need to "malloc" each gsl_vector individually and this implies that memory won't be contiguous. This hits performance because the lack of "spatial locality" in memory allocation usually increases cache misses substantially.

If you have the choice to use a completely different solution, I would implement the tensor calculation using StaticMatrix or DynamicMatrix classes in Blaze lib

Blaze

Why Blaze?

(i) The StaticMatrix or DynamicMatrix interface is much better than std::vector<gsl_vector*> or gsl_matrix

(ii) Blaze is the fastest BLAS lib available in C++. It is faster than gsl if you have available Intel MKL (remember that Intel MKL is faster than gsl BLAS). Why so? Because Blaze uses a new technique called "Smart Expression Template" . Basically, researchers in Germany showed in a series of articles paper 1 paper 2 that the "Expression Template" technique, which is the standard technique in many C++ BLAS libraries, is terrible for matrix operations (BLAS 3 operations) because compiler can't be smarter than low level code. However, "Expression Template" can be used as a smarter wrapper to low level BLAS libraries like intel MKL. So they create "Smart Expression Template" technique which is just a wrapper to your choice of low level blas lib. Their benchmarks are astonishing

benchmark

Vivian Miranda
  • 2,467
  • 1
  • 17
  • 27
  • Funny how is benchmark results shown by project X, it so happens that X comes out fastest. The Eigen folks are also quite adamant that their code is good, and Armadillo has changed / improved since the version used in the benchmark. Your mileage may vary. – Dirk Eddelbuettel Aug 23 '13 at 18:36
  • It is much faster than Armadillo (you can see the benchmark link) because as I said blaze use "Smart Expression Template". Meaning: Blaze uses templates only as a wrapper to low level c code blas libraries (like intel mkl). Armadillo tries to calculate matrix operation using loop unrolling at compile time (the original form of "expression template") and the papers I linked showed that this is far from the optimal solution (you can't beat intel mkl using just the compiler!). I think that Eigen uses a mixed implementation, but Blaze is also faster than Eigen. – Vivian Miranda Aug 23 '13 at 18:41
  • Yes. Benchmark may vary. But the idea that "Expression Template" can't beat intel MKL (or any other c low level code) is very solid. Using templates as a wrapper to optimized low level code is much better than relying on compiler to do 100% of the work. I think that all C++ BLAS are converging on that, which is the idea behind of "Smart Expression Templates". – Vivian Miranda Aug 23 '13 at 18:47
  • MKL happens to be multicore, but it is no longer the only multicore BLAS library. Armadillo calls out to BLAS as well, so can use MKL or OpenBLAS or whatever else you have installed. Eigen generally tries to its own thing rather than calling BLAS. But thanks for the pointer to Blaze, I'll keep an eye on it. – Dirk Eddelbuettel Aug 23 '13 at 18:49