[ublas] Efficiency of product of two sparse matrices

Subject: [ublas] Efficiency of product of two sparse matrices
From: George Slavov (gpslavov_at_[hidden])
Date: 2012-12-04 16:02:32

In my project I need to form the product A^T A where A is a
compressed_matrix, but it's proving to be a big bottleneck in my code,
taking on the order of 10 seconds. The result is also very sparse. The
matrix has a size of about 13000x13000 and has about 37000 nonzero
entries. The product of sparse matrices is generally not sparse, but
my matrix is structured enough that the product has roughly the same
sparsity as A.

The code I'm using is basically.

typedef ublas::compressed_matrix<T, ublas::column_major, 0,
ublas::unbounded_array<unsigned int>, ublas::unbounded_array<T> > sparse_matrix;

sparse_matrix A(13000, 13000)
sparse_matrix result(13000, 13000)
// fill A
sparse_prod(trans(A), A, result, false)

Are there any conditions on the use of sparse_prod which I'm not aware
of with regard to matrix orientations? My A is column major so the
transpose is row_major. I would have thought that would be the optimal
storage to form the product since sums range over rows of A^T and
columns of A.

Regards,
George