Re: [boost] [ublas] Weakness in large sparse matrices

From: Johan Jansson (johanjan_at_[hidden])
Date: 2003-02-17 06:58:11

On Sat, 8 Feb 2003, Joerg Walter wrote:

> Hi Johan,
>

[...]

> On 32-bit systems the upper dimensions for sparse_matrix are 65535 x 65535
> due to the internal mapping of (i, j) -> (i * size2+j) (or (i + j * size1).
> I've added a corresponding check to the debug version.
>
> For larger dimension consider to use compressed_matrix (FORTRAN CRS/CCS
> format) or coordinate_matrix (FORTRAN COO format) please. I'm also assuming
> you're checking the boost_1_29_0 release, so you'll probably have to look
> into the current CVS version for these (many changes).
>

[...]

>
> Using the experimental (means undocumented ;-( axpy_prod() functions one
> achieves the correct complexity.
>

[...]

>
> Alexei Novakov's benchmarks show an overhead of around 100% for the better
> ublas sparse matrix vector multiply implementations currently. He proposed
> an improvement of sparse matrix iterators already, which we'll  tackle
> immediately after the upcoming 1_30_0 release.
>
> Thanks for your feedback,
>
> Joerg
>
>

Hi Joerg,

Thanks for all your info. I've run the tests with a Boost from CVS (from
january 31st), compressed_matrix and axpy_prod, and the results give
roughly the same speed as our implementation, and ca. 30% better memory
efficiency. Great! The -DNDEBUG flag also seems critical, without it
performance is terrible (quadratic).

Alexei's proposed optimizations seem interesting. I tried the axpy_prod
you provided, but it didn't give any significant change. I trust your
figures however.

I will propose that we start using ublas as soon as the linear complexity
functions appear in the stable branch.

I provide our benchmark below for reference (with the timing calls, and
other dependencies stripped out):

Thanks,
  Johan

---
#include <iostream>
#include <boost/numeric/ublas/config.hpp>
#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/vector_sparse.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/matrix_sparse.hpp>
#include <boost/numeric/ublas/io.hpp>
#include <boost/numeric/ublas/operation.hpp>
#define N 1000000
using namespace boost::numeric::ublas;
namespace ublas = boost::numeric::ublas;
typedef boost::numeric::ublas::compressed_matrix<double> uBLASSparseMatrix;
typedef ublas::vector<double> uBLASVector;
int main()
{
  cout << "Benchmark for linear algebra in uBLAS" << endl;
  cout << "--------------------------------------" << endl;
  cout << "Creating a vector with " << N << " elements" << endl;
  uBLASVector x(N);
  cout << "Creating another vector with " << N << " elements" << endl;
  uBLASVector y(N);
  cout << "Creating a " << N << " x " << N << " matrix" << endl;
  uBLASSparseMatrix A(N, N);
  cout << "Assembling" << endl;
  for (int i = 0; i < N; i++) {
    A(i,i) += 2.0;
    if ( i > 0 )
      A(i, i - 1) += -1.0;
    if ( i < (N-1) )
      A(i, i + 1) += 1.0;
  }
  cout << "Multiplying: y = A*x" << endl;
  axpy_prod(A, x, y);
}
---