Thomas Willhalm wrote:

I'm new to ublas, so it's easily possible that I overlooked something obvious. The first thing I tried to calculate was the following line:

numerics::matrix<double> M(dimension,n); numerics::matrix<double> Covariance = numerics::prod(M,numerics::trans(M))/n;

where dimension is 50 and n is 8000.

Unfortunately. it turned out that the following "stupid" code actually performs better:

numerics::matrix<double> Covariance(dimension,dimension); for (int i=0; i<dimension; ++i) for (int j=0; j<dimension; ++j) { Covariance(i,j)=0; for (int k=0; k<n; ++k) Covariance(i,j) += M(i,k)*M(j,k); Covariance(i,j) /= n; }

What am I doing wrong?

I'm using gcc 2.95.3 under Linux 2.4.18 on a mobile Pentium III, if it matters.

First of all, it seems that you have fairly old version of ublas. New version is in namespace `boost::numeric::ublas'. But I don't think that this influences the performance ;o) You must define NDEBUG. Otherwise expression templates are not enabled. You can also try to compile with -O2 or -O3. For the following test program, compiled with gcc 3.2, I got: without -DNDEBUG: ublas: 9.36 loop: 5.74 with -DNDEBUG: ublas: 5.6 loop: 4.37 with -DNDEBUG -O3: ublas: 0.9 loop: 1.59 with -DNDEBUG -O3 -funroll-loops: ublas: 0.81 loop: 1.6 Test program: ============================================= #include <iostream> #include <boost/timer.hpp> #include <boost/numeric/ublas/matrix.hpp> int main() { int dimension = 50; int n = 8000; boost::numeric::ublas::matrix<double> M (dimension, n); for (int i = 0; i < dimension; ++i) for (int j = 0; j < n; ++j) M(i,j) = i+j; boost::timer t; boost::numeric::ublas::matrix<double> Covariance1 = boost::numeric::ublas::prod (M, boost::numeric::ublas::trans (M)) / n; std::cout << "ublas: " << t.elapsed() << std::endl; t.restart(); boost::numeric::ublas::matrix<double> Covariance (dimension, dimension); for (int i=0; i<dimension; ++i) for (int j=0; j<dimension; ++j) { Covariance(i,j)=0; for (int k=0; k<n; ++k) Covariance(i,j) += M(i,k)*M(j,k); Covariance(i,j) /= n; } std::cout << "loop: " << t.elapsed() << std::endl; } ==================================================== Sincerely, fres