You're comparing statically-defined arrays which live on the stack:
to heap allocated arrays:
mat a1= randi<mat>(100, 100, distr_param(1, 10));
mat b1= randi<mat>(100, 100, distr_param(1, 10));
Of course, stack-allocated variables are ...
Nested used of OpenMP is allowed. However, by default the second level will only get a single thread. If that's ok with you, you're good now. If you want the embedded level to be parallel too, set the environment variable OMP_NESTED=true or a function call omp_set_nested(1).
The answer to your question is benchmarking. Both Armadillo and Eigen provide some benchmarks in their documentation. The only problem is that they don't compare to the same libraries, but you can still get an idea. Maybe there are more comprehensive benchmarks elsewhere...
In any case, the overhead of wrapping LAPACK in a C++ library will be much smaller ...
So I invented my own solution. Critique, suggestions are highly appreciated.
I learned that elements are stored in a matrix as column-major. So I assumed that using std::copy() with an iterator over a single column will just give consecutive elements, hence better speed.
Here's my solution:
template <typename T>