Hi!

I am new to Eigen2. What concerns me is the interoperability with other libs, such as uBlas.

In my project I make use of sparse and dense linalg objects. I am planning to start porting parts of the project to eigen2 while keeping uBlas at those places where it provides features Eigen2 does not support yet. Therefore, I often efficiently need to apply linalg operations onto mixed eigen2/ublas objects.

uBlas solved bindings to other linalg packages by providing a double* data array which is filled according to the predefined storage format (eg. sparse row_major, dense column_major) . Frankly, I have not seen in the API doc how I can swap the contents of an Eigen2 sparse_matrix with an uBlas compact-matrix or dense vectors. (eg. read data from BLAS data/ export data to BLAS) I also miss a function where I can efficiently fill a sparse matrix (for example, if I already know the sort order of the input data and want to avoid the sort() after an item was inserted).

In the case that the features I am asking for are not yet available, is there a special version of Eigen2 when I can expect them?

Thanks

Seb wrote:In my project I make use of sparse and dense linalg objects. I am planning to start porting parts of the project to eigen2 while keeping uBlas at those places where it provides features Eigen2 does not support yet. Therefore, I often efficiently need to apply linalg operations onto mixed eigen2/ublas objects.

Sounds reasonable.

uBlas solved bindings to other linalg packages by providing a double* data array which is filled according to the predefined storage format (eg. sparse row_major, dense column_major) . Frankly, I have not seen in the API doc how I can swap the contents of an Eigen2 sparse_matrix with an uBlas compact-matrix or dense vectors. (eg. read data from BLAS data/ export data to BLAS) I also miss a function where I can efficiently fill a sparse matrix (for example, if I already know the sort order of the input data and want to avoid the sort() after an item was inserted).

For dense matrices, Eigen lets you access their array of data thanks to the data() method (see in Matrix). Eigen lets you map any existing array as a Eigen expression: see class Map and the static methods Map() and MapAligned() in Matrix. That's very useful in cases like yours where you want to mix Eigen and non-Eigen stuff.

About the sparse stuff, I really don't know, Gael is the only reference here. I hope he'll help you here, otherwise you can find his e-mail in the source code (like Eigen/src/Core/MatrixBase.h).

At least I know that there are startFill() and endFill() methods in SparseMatrix that sound like what you're looking for.

Hi Seb,

I don't know very well ublas, but if I remember well I know they provide various fancy ways to store matrix data. So for dense matrix, assuming you use ublas matrix types for which data are sequentially stored in memory (row-major or col-major, Eigen can handle both), then the Eigen::Matrix::Map and Eigen::Matrix::data()
functions are what you are looking for. Some examples:

Code: Select all

Eigen::MatrixXf eigenmat = Eigen::MatrixXf::Map(ublasmat.data(),
ublasmat.rows, ublasmat.cols);

or you can also directly use a Map object that has the advantage to
avoid the copy of the data:

Code: Select all

Eigen::Map mappedObject(ublasmat.data(),
ublasmat.rows, ublasmat.cols);
mappedObject.lu().solve(...);

For the other way round, mat.data() gives you a pointer to the data, and I guess ublas provide the feature to import such data.

Now, if you want, you can go a bit further by extending the MatrixBase and Matrix classes. See: http://eigen.tuxfamily.org/api/Customiz ... MatrixBase. Those explanation are for MatrixBase, but you can do the same for Matrix (the name of the macro to define is EIGEN_MATRIX_PLUGIN). For instance you can overload operator= for ublas objects, add a convenient toUblas() function, etc... If you need help to do that, feel free to ask.


About sparse matrix. First of all, I must say that the sparse module is not ready yet, even though there are already a lot of stuff. For efficiency and compatibility reasons, so far we support only one kind of sparse matrix storage: "compressed row/col storage" which is used by many other libs. In Eigen, all linear algebra operations on sparse matrix always return a sparse matrix with sorted coefficients. To fill a sparse matrix, the best way is to fill it in fixed order such that there is no need to sort the data, e.g.:

Code: Select all

SparseMatrix mat(rows,cols); // the default is col-major
mat.startFill(  );
for (int j=0; j )
        mat.fill(i,j) = ;
mat.endFill();

If you cannot fill your matrix in such a way, we have the fillrand() function which is more flexible: it only requires the matrix is filled by increasing outer index, eg:

Code: Select all

SparseMatrix mat(rows,cols); // the default is col-major
mat.startFill(  );
for (int j=0; j;
  }
mat.endFill();

This method is still very fast when you have only a few nonzeros per column/row, let's say up to 32. If your write accesses are completely random, then we provide a RandomSetter which temporarily store the data using a hash-map (the best is to use googles's dense hash map class, we provide a support for it).

If you plane to use the sparse module, I highly recommend you to use svn trunk (and generate the doc from it). I'll try to extend the documentation as soon as I can.

Also, I'd be very interested to know what are the ublas features that Eigen does not provide yet. This might help me a lot to know on which part of the sparse module I should focus.

For your information, our sparse module already provide a large set of dense solvers with support for a couple of 3rd party libs, see: http://eigen.tuxfamily.org/index.php?ti ... rent_state

Thanks to both of you for you fast and well reply. I wish you a happy new year at this point!

The direct application of uBlas algorithms onto C-Arrays is not that easy as I expected. One has to deal with adaptor classes which are first transformed into a vector class (using an adaptor as storage format) and then convert this into the final vector.

For example, if one wants to access the Eigen2 object A via a function

Code: Select all

void funct(const ublas::vector v);

one has to call

Code: Select all

funct( ublas::vector > (6,ublas::shallow_array_adaptor(6,A.data())) );

which does anything but not fast. Maybe it is possible in the future to provide some kind of binding such that an Eigen2 object may be converted to ublas directly?

I started porting and I am impressed of the performance boost. Regarding the features, uBlas still offers in the sparse module:

1. generic types: Is it possible to create a sparse object of Eigen::SparseMatrix ?
2. I can not multiply A=B*C where A and B are sparse matrices. Currently, A must be dense.

For the future plans I'd like to see the following features in Eigen2 (I often deal with equation systems arising from variational problems that are symmetric):

1. special products:
   v'Av, B'AB where A is a symmetric dense/sparse marix, B a sparse/dense matrix and v a vector.
2. symmetric matrix-vector products:
   Av where A is symmetric dense/sparse matrix. These 2 wishes arise when dealing with optimization problems and constraints.
3. sparse matrix types with columns/rows being known at compile time (may this speed up something as well? It should reduce the size of small sparse objects being allocated and thus reduce computation time)
4. sparse matrix types which are specialized on very small and very large dimensions. I have seen that currently int is used for indices. This may not satisfy my need in case of very large matrices (at least for the number of nonzeros). But I also agree that chosing size_t is not recommended in case of very small matrices when (losely spoken) the allocation of internal variables probably takes more memory than the data storage.

For the special products I currently compute
D=BAB' by C=BA and D=C*B' via B.transpose(); B=sparse. Is this the most efficient way at the moment? I wonder what transpose() is doing and if it effects the multiplication...

Thanks for listening.

Seb wrote:Thanks to both of you for you fast and well reply. I wish you a happy new year at this point!

The direct application of uBlas algorithms onto C-Arrays is not that easy as I expected. One has to deal with adaptor classes which are first transformed into a vector class (using an adaptor as storage format) and then convert this into the final vector.

For example, if one wants to access the Eigen2 object A via a function

Code: Select all

void funct(const ublas::vector v);

one has to call

Code: Select all

funct( ublas::vector > (6,ublas::shallow_array_adaptor(6,A.data())) );

which does anything but not fast. Maybe it is possible in the future to provide some kind of binding such that an Eigen2 object may be converted to ublas directly?
[/code]

Well, if it appears that ublas matrix/vector store their attributes like (1 or 2 int for the dimension(s), and 1 pointer for the data) Eigen's Matrix does, then you can easily add a asUblas() function as a Matrix addons which would simply do a reinterpret_cast towards the respective ublas type. Some meta programming will be needed to retrieve the equivalent ublas type.

I started porting and I am impressed of the performance boost. Regarding the features, uBlas still offers in the sparse module:

1. generic types: Is it possible to create a sparse object of Eigen::SparseMatrix ?
2. I can not multiply A=B*C where A and B are sparse matrices. Currently, A must be dense.

For the first one, yes, Eigen::SparseMatrix should work but you have to make sure that myObjectType support common operators (=,+,-,*,etc.). You will also need to overload a couple of other functions such as ei_abs(myObjectType). See http://eigen.tuxfamily.org/api/Customiz ... ScalarType.

Actually, A=B*C *must* work with A, B and C sparse matrices. We have unit tests doing exactly that.

[*] special products:
v'Av, B'AB where A is a symmetric dense/sparse marix, B a sparse/dense matrix and v a vector.

For the case of vectors, I don't think there is much to do for that special cases. For B'AB, yeah, maybe, do you know if some libraries have optimized routines for them, that would be interesting to bench.

[*] symmetric matrix-vector products:
Av where A is symmetric dense/sparse matrix. These 2 wishes arise when dealing with optimization problems and constraints.

There is not much to do. Because of SSE optimization, for dense matrices, it is much faster to consider to store all coeffs of A and treat A as a standard matrix. And, well if A is sparse, I still don't see what could be optimized (excepted the memory by storing only the upper/lower part as the price of more expensive computation for the product).

[*] sparse matrix types with columns/rows being known at compile time (may this speed up something as well? It should reduce the size of small sparse objects being allocated and thus reduce computation time)

Hm, small sparse objects does not make sense: use a plain dense matrix for that, that will always be much much faster, and a small sparse object is likely to require more memory than a plain dense object.

[*] sparse matrix types which are specialized on very small and very large dimensions. I have seen that currently int is used for indices. This may not satisfy my need in case of very large matrices (at least for the number of nonzeros). But I also agree that chosing size_t is not recommended in case of very small matrices when (losely spoken) the allocation of internal variables probably takes more memory than the data storage.

Yeah, I guess I could template the type of the indices. But do you really think you might have matrices with more than 2^32 non zero coefficients ? That's really huge. For instance, with 2^32 nonzeros of double and using int64 for the indices would need 64GB of memory. I think that for such problems you need big clusters with highly parallelized libraries.

For the special products I currently compute
D=BAB' by C=BA and D=C*B' via B.transpose(); B=sparse. Is this the most efficient way at the moment? I wonder what transpose() is doing and if it effects the multiplication...

I'd recommend to let eigen do its job with the temporaries and simply write:
C = B.transpose() * A * B;

.transpose() has zero direct cost, but indirectly it changes the way the product is performed. Sometime, for sparse matrices we have to copy a term to another storage order such that the product can be done efficiently. Fortunately, Eigen's transposition of sparse matrix is very fast.
The slow operations are:
row-major = col-major * col-major
and
col-major = row-major * row-major

I am sorry for repeatingly bother you with stupid questions - I am not aware of Eigens internals...

For the first one, yes, Eigen::SparseMatrix should work but you have to make sure that myObjectType support common operators (=,+,-,*,etc.). You will also need to overload a couple of other functions such as ei_abs(myObjectType). See http://eigen.tuxfamily.org/api/Customiz ... ScalarType.

Actually, A=B*C *must* work with A, B and C sparse matrices. We have unit tests doing exactly that.

It did not work with the last official beta, but good to hear that it is already there! (I have some problems with the current svn - see below)

[*] special products:
v'Av, B'AB where A is a symmetric dense/sparse marix, B a sparse/dense matrix and v a vector.

For the case of vectors, I don't think there is much to do for that special cases. For B'AB, yeah, maybe, do you know if some libraries have optimized routines for them, that would be interesting to bench.

The case of vectors is, indeed, very common. The case of dense matrices appears when transformations/reductions are used in dynamical simulations. Therein one linearly transforms the variable space of some energy function using dense trafo matrices B where A is the 2nd energy derivative with respect to the original variables.
The case of B being a sparse matrix appears in most finite element codes only once per runtime, namely by applying linear constraints via static condensation. Unfortunately I do not know of any efficient code doing that. I read some items on uBlas' mailing list where it was suggested to perform the operation vector wise. The problem is more or less to deal with the transpose of B and with the temporary which may be rather dense, but must be stored in sparse form in order to save memory. Off course, I could live with 2 subsequent products as long as they are sufficient fast.

[*] symmetric matrix-vector products:
Av where A is symmetric dense/sparse matrix. These 2 wishes arise when dealing with optimization problems and constraints.

There is not much to do. Because of SSE optimization, for dense matrices, it is much faster to consider to store all coeffs of A and treat A as a standard matrix. And, well if A is sparse, I still don't see what could be optimized (excepted the memory by storing only the upper/lower part as the price of more expensive computation for the product).

Thanks for explaining the dense case. For sparse: The memory is exactly the problem when considering sparse FEM matrices with rather many nonzeros. Why should one keep more than required?

[*] sparse matrix types with columns/rows being known at compile time (may this speed up something as well? It should reduce the size of small sparse objects being allocated and thus reduce computation time)

Hm, small sparse objects does not make sense: use a plain dense matrix for that, that will always be much much faster, and a small sparse object is likely to require more memory than a plain dense object.

On the level of material laws one often deals with objects of dimension 3,6 or 9. The corresponding matrices are often symmetric and very sparse (fill ratio 1:4...1:6). Coming from uBlas I figered out that using sparse objects (eg. 6x9) multiplied with 6x6 dense matrices was significantly faster than using dense counterparts (in that special case). If time is left I have to test this again after porting to Eigen.

Yeah, I guess I could template the type of the indices. But do you really think you might have matrices with more than 2^32 non zero coefficients ? That's really huge. For instance, with 2^32 nonzeros of double and using int64 for the indices would need 64GB of memory. I think that for such problems you need big clusters with highly parallelized libraries.

Oops - I still can't count magnitudes... Though I sometimes touch the 32G limit, it still is representable...

I'd recommend to let eigen do its job with the temporaries and simply write:
C = B.transpose() * A * B;

Thanks.

The slow operations are:
row-major = col-major * col-major
and
col-major = row-major * row-major

I did not yet get into when Egen2 is using row-major and when not. Is it chosen automatically?

Hi Gael,
I have problems filling a sparse matrix in current SVN?
Following your advice above, I copy the contents of an uBlas sparse matrix 'stiffness' to an Eigen2 object 'eigen_k_full' and a triangular matrix 'eigen_k_tri':

Code: Select all

   eigen_k_full.startFill(stiffness.nnz());
   eigen_k_upper.startFill(stiffness.nnz()/2);
   for (sfem::dofs::DofHandler::Type_sparse_matrix_assemble::iterator1 i1 = stiffness.begin1(); i1 != stiffness.end1(); ++i1)
   {
      int i,j;
      for (sfem::dofs::DofHandler::Type_sparse_matrix_assemble::iterator2 i2 = i1.begin(); i2 != i1.end(); ++i2)
      {
         double s=*i2;
         if (s)
         {
            i=i2.index1();
            j=i2.index2();
            eigen_k_full.fillrand(i,j) = s;
            if (i<j) eigen_k_upper.fillrand(i,j) = s;
            if (i==j) eigen_k_diag(i) = s;
         }
      }
   }
   eigen_k_full.endFill();
   eigen_k_upper.endFill();


The dimensions are 169983x169983. Then I counted the nonzeros:

Code: Select all

A = 6140203 (nnz uBlas)
B = 118529 (nnz Eigen2), should be A
C = 1319 (nnz Eigen2_triangular), should be around A/2


Moreover, I tried a Matrix-vector multiplication for the cases A and B:

Code: Select all

UBLAS 00:00:00.081642
EIGEN 00:02:45.202432


Seb wrote:I am sorry for repeatingly bother you with stupid questions - I am not aware of Eigens internals...

For the first one, yes, Eigen::SparseMatrix should work but you have to make sure that myObjectType support common operators (=,+,-,*,etc.). You will also need to overload a couple of other functions such as ei_abs(myObjectType). See http://eigen.tuxfamily.org/api/Customiz ... ScalarType.

Actually, A=B*C *must* work with A, B and C sparse matrices. We have unit tests doing exactly that.

It did not work with the last official beta, but good to hear that it is already there! (I have some problems with the current svn - see below)

hm, that's weird, that should really work. btw, what do you mean by does not work ? compilation issue, runtime error, wrong result ?

[*] special products:
v'Av, B'AB where A is a symmetric dense/sparse marix, B a sparse/dense matrix and v a vector.

For the case of vectors, I don't think there is much to do for that special cases. For B'AB, yeah, maybe, do you know if some libraries have optimized routines for them, that would be interesting to bench.

The case of vectors is, indeed, very common. The case of dense matrices appears when transformations/reductions are used in dynamical simulations. Therein one linearly transforms the variable space of some energy function using dense trafo matrices B where A is the 2nd energy derivative with respect to the original variables.
The case of B being a sparse matrix appears in most finite element codes only once per runtime, namely by applying linear constraints via static condensation. Unfortunately I do not know of any efficient code doing that. I read some items on uBlas' mailing list where it was suggested to perform the operation vector wise. The problem is more or less to deal with the transpose of B and with the temporary which may be rather dense, but must be stored in sparse form in order to save memory. Off course, I could live with 2 subsequent products as long as they are sufficient fast.

yes, I know such products are quite common, but frankly I think that doing 2 products should be fast enough, and I'd expect much gain from a specialized implementation... What would be very interesting it is know the respective nnz of B'A and B'AB. If in practice it appears that B'AB is more sparse than the temporary B'A (or AB) then that would be a good argument in favor of specialized implementation.

[*] symmetric matrix-vector products:
Av where A is symmetric dense/sparse matrix. These 2 wishes arise when dealing with optimization problems and constraints.

There is not much to do. Because of SSE optimization, for dense matrices, it is much faster to consider to store all coeffs of A and treat A as a standard matrix. And, well if A is sparse, I still don't see what could be optimized (excepted the memory by storing only the upper/lower part as the price of more expensive computation for the product).

Thanks for explaining the dense case. For sparse: The memory is exactly the problem when considering sparse FEM matrices with rather many nonzeros. Why should one keep more than required?

ok, note that the memory overhead in lower than a factor 2, so not that much. But I've though a bit about that (well, 2 secs) and actually even in the sparse case we can implement something very efficient (assuming the result of Av is quite dense). So be sure this will be done in Eigen.

[*] sparse matrix types with columns/rows being known at compile time (may this speed up something as well? It should reduce the size of small sparse objects being allocated and thus reduce computation time)

Hm, small sparse objects does not make sense: use a plain dense matrix for that, that will always be much much faster, and a small sparse object is likely to require more memory than a plain dense object.

On the level of material laws one often deals with objects of dimension 3,6 or 9. The corresponding matrices are often symmetric and very sparse (fill ratio 1:4...1:6). Coming from uBlas I figered out that using sparse objects (eg. 6x9) multiplied with 6x6 dense matrices was significantly faster than using dense counterparts (in that special case). If time is left I have to test this again after porting to Eigen.

That would be very cool if you could do some benchmark since I think Eigen's dense path is much faster than ublas one.

The slow operations are:
row-major = col-major * col-major
and
col-major = row-major * row-major

I did not yet get into when Egen2 is using row-major and when not. Is it chosen automatically?

By default it is col-major, if you want a row major matrix :

Code: Select all

SparseMatrix mat;


Seb wrote:Hi Gael,
I have problems filling a sparse matrix in current SVN?
Following your advice above, I copy the contents of an uBlas sparse matrix 'stiffness' to an Eigen2 object 'eigen_k_full' and a triangular matrix 'eigen_k_tri':

Code: Select all

   eigen_k_full.startFill(stiffness.nnz());
   eigen_k_upper.startFill(stiffness.nnz()/2);
   for (sfem::dofs::DofHandler::Type_sparse_matrix_assemble::iterator1 i1 = stiffness.begin1(); i1 != stiffness.end1(); ++i1)
   {
      int i,j;
      for (sfem::dofs::DofHandler::Type_sparse_matrix_assemble::iterator2 i2 = i1.begin(); i2 != i1.end(); ++i2)
      {
         double s=*i2;
         if (s)
         {
            i=i2.index1();
            j=i2.index2();
            eigen_k_full.fillrand(i,j) = s;
            if (i<j) eigen_k_upper.fillrand(i,j) = s;
            if (i==j) eigen_k_diag(i) = s;
         }
      }
   }
   eigen_k_full.endFill();
   eigen_k_upper.endFill();


The dimensions are 169983x169983. Then I counted the nonzeros:

Code: Select all

A = 6140203 (nnz uBlas)
B = 118529 (nnz Eigen2), should be A
C = 1319 (nnz Eigen2_triangular), should be around A/2


Moreover, I tried a Matrix-vector multiplication for the cases A and B:

Code: Select all

UBLAS 00:00:00.081642
EIGEN 00:02:45.202432


hm.. note that if eigen_k_full and eigen_k_upper are col-major (that is the default), then you must make sure that j always increase. For instance if at some point j=10, then cannot add new non zero into the column of index c with c<j. If you compile in debug mode, then you should get an assert for such errors. If you cannot set the matrix with such a coherence then you can use the RandomSetter wrapper.

You were right. Thanks for the hint to RandomSetter.

What still bothers me is the bad performance of the matrix-vector product. It got worse when using the properly filled matrix. (7:50 minutes compared with 0.1 seconds using uBlas) -> factor 4500 slower!!

Another thing is the memory usage:
Filling the two matrices simultaneosuly (upper triangular and full version) requires about 600 MB RAM. How is that possible?

For the product I used (NDEBUG mode):

Code: Select all

// UBLAS:
typedef ublas::vector > vector_type;
typedef ublas::generalized_vector_of_vector > > sparse_assemble_matrix;
sparse_assemble_matrix stiffness;
vector_type vec01 = ublas::prod( stiffness, vector_ublas );
// EIGEN2:
Eigen::VectorXd eigen_vec = Eigen::VectorXd::Zero( numDofs );
Eigen::VectorXd vec02 = eigen_k_full*eigen_vec;


I could save the input matrix to a text file in order to send it to you.

EDIT: The memory measure was losely done using 'top'. It may cover assembly and product as well.

You recommended to use dense matrices disrespecting their known dense or sparse structure if the dimensions are very small. I can now confirm your suggestion.

Eigen2 sparse matrices are - indeed - very slow at the moment. After porting most uBlas stuff to Eigen2 I compared some times in the small-object module (where 6x9 sparse matrices are used). When using these matrices where dense_marix(6,6) * sparse_matrix the computing times increase by a factor of at least 4 for the overall algorithm. When using Eigen2 where the product sparse_matrix(9,6)*vector(6,1) appears, the overall algorithm still improves performance by a factor of almost 2. I regard this boost of speed to Eigen2 dense operations where uBlas is known to be quite slow.

I'd suggest to benchmark Eigen2 against uBlas' sparse package. The differences in computing time and memory usage are too obvious right now...

Seb wrote:You recommended to use dense matrices disrespecting their known dense or sparse structure if the dimensions are very small. I can now confirm your suggestion.

Eigen2 sparse matrices are - indeed - very slow at the moment. After porting most uBlas stuff to Eigen2 I compared some times in the small-object module (where 6x9 sparse matrices are used). When using these matrices where dense_marix(6,6) * sparse_matrix the computing times increase by a factor of at least 4 for the overall algorithm. When using Eigen2 where the product sparse_matrix(9,6)*vector(6,1) appears, the overall algorithm still improves performance by a factor of almost 2. I regard this boost of speed to Eigen2 dense operations where uBlas is known to be quite slow.

Wow. Your matrix are 6x9, that's really really small!!

Is this size 6x9 fixed at compile-time? (From what you say, it sounds so).

If yes, you should really used a fixed-size Eigen matrix type here, it will improve your performance considerably.

So instead of using "MatrixXd my_matrix(6,9);" you could try "Matrix my_matrix;"

Likewise many Eigen methods have fixed-size variants, be sure to use them whenever applicable, the speed gain is really large.

fixed-size matrices are fully compatible with dynamic-size ones, with Eigen.

Side note: if on the other (as what you say suggests) you have 9x6 matrices (as opposed to 6x9), and you use "double" numbers, it may help (for better vectorization) to use row-major storage, as then the rows will then be 6*sizeof(double)==48 bytes which is a multiple of 16 bytes which means that vectorization will be possible in this direction.

Seb wrote:What still bothers me is the bad performance of the matrix-vector product. It got worse when using the properly filled matrix. (7:50 minutes compared with 0.1 seconds using uBlas) -> factor 4500 slower!!

it is simply because the Sparse module is far to be ready yet. So we have sparse * sparse, but the sparse * dense is not properly handled yet. So what happens is that it tales the dense path, so it performs a dense mat * vec product where for each access to the matrix it performs a search of the coefficient...

Now handling sparse_mat * dense_vector is trivial because internally, sparse * sparse is able to use a dense vector is the result appears to be quite dense. What is missing is a little of meta-programming around....

seb:

sparse matrix * vector is now optimized. I get better perf than MTL4 and GMM++. I hope that works for you too.

also, can you tell us exactly which operations are slow. Because I consistently ran some benchmark of the sparse module against GMM++, MTL4 and CSparse, and Eigen is always faster for the operations I tried. So, this is is what I benched so far:

* linear stuffs like add, sub, etc.
* sparse = sparse * sparse, sparse = sparse * sparse.transpose(), sparse = sparse.transpose() * sparse, sparse = sparse.transpose() * sparse.transpose()
* dense_vec = sparse * dense_vec, dense_vec = sparse.transpose() * dense_vec
* dense_vector = sparse.triangularSolve(dense_vector)

all these benchmarks are in eigen2/bench.

So, if you find an operation which is really slower than ublas, please tell us.

Gael,

thanks for this update!! Sounds very nice! Thanks for the effort! Frankly, I won't be able to test it immediately.

I was thinking about a comment of yours on sparse*sparse products: When developing Finite Element solvers, one often has to deal with constraint equations. As an example, a 2nd derivative of the energy function with say dimension=1e6 and many nonzeros is paired with a set of constraint equations whose gradients are very sparse, say each constraints depends on 5..50 variables, then it would be, indeed, very interesting to have at least the choice between two sparse=sparse*sparse products that are either optimized using dense vectors or sparse arrays... The default method ('*') may be tuned to rather dense vectors, which sounds reasonable for most cases.

To come back onto the symmetric sparse matrices: Right now I implemented the assemble process of a matrix A by assembling B with A=B+B' using uBlas. The product f=A*v using simple uBlas routines Bv+B'v is considerably slower than iterating over the nonzeros and doing for each iterate 'it' f(it.i)+=it.value*v(it.j); f(it.j)+=it.value*v(it.i) which is not slower than the product of the full matrix. The memory usage decreased significantly, though.