[SciPy-user] Sparse factorisation

[Gael Varoquaux]

On Thu, Apr 30, 2009 at 11:47:38PM -0400, Nathan Bell wrote:
> On Thu, Apr 30, 2009 at 11:54 AM, Gael Varoquaux
> <gael.varoquaux at normalesup.org> wrote:

> > I have dug a bit further. It seems that all that needs to be done is
> > expose the L, U, perm_c and perm_r attributes of the SciPyLUObject (AKA
> > factored_lu object) returned by splu. Of course this easier said than
> > done, as exposing these objects requires creating scipy sparse matrices
> > from the inner SuperLU representation of these objects.

> > I'd really appreciate if someone (most probably Nathan) could give me a
> > hand with this. I realise that I am asking people to do my work, and I
> > know exactly what my reaction is when someone comes around to me with
> > this request (ie, not happy), but I am not sure I have to time required
> > to learn the library and the tools to do this myself, and would hate to
> > have to find an ugly workaround (this does sound like a bad excuse).

> I hate to disappoint, but I don't have enough free time right now.  If
> you or someone else wants to give it a try I can probably provide some
> support.

Fair enough. I can understand that.

In the mean time, I figured out that I only needed to expose thee
permutation vectors from the factored_lu object. These permutation vector
are not sparse objects, so it is very easy to expose them to numpy (just
a question of PyArray_SimpleNewFromData and using the 'base' trick (
http://blog.enthought.com/?p=62 ) to do the reference counting. 

This is a bit ugly, because it does not expose all the information.
Another issue is that as long as there are references on the arrays
created as views from the permutation vectors, the whole factored_lu
object (much heavier in memory) is not garbaged collecter.

If we can live with these limitations, I think I can contribute a clean
patch. It does not solve completely my problem, but still does one step
in the right direction.

Gaël