Re: [eigen] SparseLU generates dense representation somewere

[Gael Guennebaud <gael.guennebaud@xxxxxxxxx>]

Handling true sparse right-hand sides in sparse solvers is much more difficult than it might look, and most cases, when solving Ax=b with both A and b sparse, x is likely to be dense. It would be interesting to have some statistics on your real-world example, that is, what are the sizes of A and b, and what are their number of non-zeros?

Moreover what's the application? I mean do you really have to compute x explicitly? For instance, in hybrid sparse solvers, the Schur complement C can be defined as the solution of a sparse problem A C = B with A and B sparse, but nobody try to compute C explicitly, but instead every time we apply C to some vectors, we apply B and then use an iterative sparse solver for applying A^-1. (the actual _expression_ of the Schur complement if more complicated, but that's the idea).

cheers,

Gael

On Sun, Jan 19, 2014 at 3:32 AM, Wenzel Jakob <wenzel@xxxxxxxxxxxxxxxxxxxx> wrote:

Thank you, that fixed the crash. But this way it's certainly too slow to be
usable, since the dense matrix is still generated slice by slice :)
(5 minutes, still running).

Do you think that it would bepossible to hack the defaultEvalTo() function
to use a SparseVector with NbColsAtOnce==1, or is it more involved than
that? At least in principle I think it would be useful if a fully sparse
solve() function is provided as an alternative option.

Thank you,
Wenzel

On 1/19/14, 3:15 AM, Christoph Hertzberg wrote:
> On 18.01.2014 22:19, Wenzel Jakob wrote:
>> So it seems like for some reason it is allocating a dense
>> representation. What am I doing wrong?
>
> It appears that the sparse solver currently takes the rhs 4 columns at a
> time, converts it to a dense matrix, solves the dense system and copies a
> sparse version of the result to the result matrix.
>
> The temporary matrix had as many columns as the result matrix, which I
> fixed here: https://bitbucket.org/eigen/eigen/commits/b75890d
> That resolves the issue of bad::alloc, but still is extremely slow for your
> simple example.
>
> I assume, if all the systems are very sparse, an optimized version which
> only uses sparse matrices/vectors could be implemented. However, that will
> likely turn out slower for matrices with higher fill grades.
>
> Christoph
>
>