Hi,

I need help with using sparse matrices in Eigen. I am writing an application for simple FEM modeling (demo on YouTube).

Currently I'm using MatrixXf class and its lu().solve(...) methods. Unfortunately this method is suitable only for small number of nodes (in my case up to ~100 on 2.66 GHz machine).

I've read tutorial on sparse matrices. Alas instructions on how to solve equations using them is still to be written.

I suppose the most suitable method for my application is SuperLU one. But using it (slu.solve(b, &V);) cause firing an assertion that says this method is not implemented (but 'test/sparse_solvers.cpp' file uses it).

How to solve equations with sparse matrices?

Side note:
In one of my projects (unfortunately w/o eigen) TausCS did a brilliant job solving FE Problems.
I would love to see a comparison of all the diffenernt sparse matrix solver, but i'm too lazy to do it myself...

hm this should really work. So first of all make sure you are using the devel branch which is highly recommended for the Sparse module. Then here is code snippet which should do the job:

Code: Select all

SparseMatrix<double> m;

// fill m

SparseLU<EigenSparseMatrix,SuperLU> slu(m, Eigen::ColApproxMinimumDegree);
bool ok = slu.solve(b,&x);

where b and x are dense vectors or matrices. Currently they cannot be sparse objects because SuperLU does not support that.

FMD, thanks for info about TausCS. I'd use it as my last resort.

ggael, Thank you for your reply. There must be something I'm doing wrong.

I got the devel branch from this link, but there is a clearly assertion "not implemented yet". So do in bitbucket "tip" tag SparseLU.h.

Code: Select all

template<typename MatrixType, int Backend>
void SparseLU<MatrixType,Backend>::compute(const MatrixType& a)
{
  ei_assert(false && "not implemented yet");
}

yes there is no default implementation for the SparseLU class and you must use one of the backends: currently SuperLU or UmfPack:

Note the second template parameter in the declaration of the sparse LU solver:

SparseLU<SparseMatrix<double>,SuperLU> slu(m);


Also note that we also support Taucs via the SparseLLT solver:

SparseLLT<SparseMatrix<double>,Taucs> llt(m);

but this is only for symmetric matrices.

Thank you very much, ggael.

I will look at these backends

I need a sparse cholesky factorizer that can quickly update an existing factorization when you change a small number of elements of the matrix being factorized. For my application, I don't care if the factorization is LLT or LDLT. Does Eigen support such factorization-updating with any of its backends? Or does LDL::compute() and LDLT::compute always recompute from scratch?

Of all the backends, I'm only familiar with CHOLMOD, but the CHOLMOD function for updating a factorization is "cholmod_sparse_updown", and I didn't find that called anywhere in the Eigen 2.0.9 source tree.

Finally, if such operations aren't supported yet, are there plans on introducing them in the near future? I need to decide soon whether to use CHOLMOD or Eigen as the solver for my project, and would love to use Eigen.

-- Matt

Indeed, currently that's not possible, however instead of directly using Cholmod I'd like to suggest you 2 options:

1 - Use Eigen's SparseMatrix with all its advantages and when you need special Cholmod features, you can use the CholmodSupport features to see a Eigen::SparseMatrix as a Cholmod matrix.

2 - Extend the Cholmod backend for SparseLLT to allow fast updates, and then I'd please to add it upstream. I don't know how does the cholmod_sparse_updown work but I guess that's only a matter of a very few lines of code to extend SparseLLT for it. If you go for that option, have a look at the src/Sparse/CholmodSupport.h file. The is a specialization of the SparseLLT class for Cholmod. The implementation of the compute() method will give you an idea how the cholmod functions are called on Eigen's sparse matrices. Then you could add an update() function doing exactly what you want. And also feel free to improve this backend as you like. E.g., by default the reordering is disabled while it is probably better to make it uses AMD by default...

Oh, in case you wonder how to use the SparseLLT class with Cholmod, have a look at the test/sparse_solvers.cpp file. There is one example per currently supported solver.

Mes wrote:I will look at these backends

Small update to the topic. I have checked some backends but none of them seemed easy to use. Being completely helpless I wrote my Vector and MatrixSparse (using CRS) classes and wrote pretty straightforward Conjugate Method solver, which is nicely presented in "Painless Conjugate Method" by Shewchuk.

Now I can easily operate with thousands of nodes in my application (analysis takes about 0.5 s, which is pretty cool).

http://dl.dropbox.com/u/460712/half-cell.png

I am very interested by the eigen project and I would like to use it in some industrial software. I am ok with using Cholmod or taucs or umfpack or whatever to direclty inverse sparse matrices, but any serious sparse direct solver rely heavily on blas 2 and 3 to get high performances (for multifrontal or supernodal methods).
So in one way or another one as to resort to a blas library instead of eigen !

Could there be an interface in eigen to blas that would permit to link some direct solver with eigen instead of a blas lib (since sounds weird since it mix C and C++ )? No chance that in a near future sparse solvers libraries rely on eigen instead of blas, so...

To me, eigen is better than any blas library (tell me if I'm wrong) because of its portability. For the moment I don't care about taking advantage of parallel hardware, but I'm sure a day eigen will do this also.

Ah, a BLAS lib implemented using Eigen? Gael has started that, see the blas/ directory in the development branch. There is much of levels 1 and 3, with level 2 missing. But this isn't our top priority at the moment as we are focusing on readying Eigen 3.0.

I posted to the wrong thread. Please ignore it