[SciPy-dev] Sparse BLAS in scipy

[Anand Patil]

On Nov 28, 2007 4:06 PM, Nathan Bell <wnbell at gmail.com> wrote:

> On Nov 27, 2007 4:35 PM, Anand Patil <anand.prabhakar.patil at gmail.com>
> wrote:
> > That and the analogous forward multiplication with a triangular matrix
> are
> > the two things I really need right now, yeah. The backend would be
> simple,
> > yet painful; upper/lower, transpose and side arguments, in addition to
> > accounting for all the possible combinations of matrix formats, and the
> > strides when 'B' is dense, would make for a huge number of bug
> > opportunities. That's what made me think wrapping a library would be
> less
> > work than writing even 'two' routines from scratch.
>
> I guess I don't fully understand what you need here.  What I imagined
> was two methods for forward/backward solves with lower/upper
> triangular matrices in CSR format.
>

I was thinking the triangular matrices could be any format, and the 'b'
matrix/vector could be any format, in which case the multiplicity of
functions rivals that of the ordinary multiplications involving sparse
matrices.

BTW, for such matrices, isn't this is equivalent to Gauss-Seidel[1]
> sweep in the appropriate direction.  If so, then we might just move
> something like [2] to scipy.linsolve.
>

I can't see any difference... it would be nice to have the L3 version too,
though, and in that case b could be in any format.

Also, are you sure the sparse LU solvers wouldn't do this for you anyway?
>

Nope, will have to look more closely at those.


> About a year ago I implemented sparsetools from scratch.  I looked at
> the same packages you listed before and came to the conclusion that
> the amount of code needed was rather small and that by writing the
> library to fit well with SciPy was worth the effort.
>
> Ultimately, I think making the backend code amenable to SWIG and
> scipy.sparse saved more effort than starting with an existing library.


OK, if you've done the research already and come to that conclusion I should
just do it! Regardless of the backend (LU, gauss-seidel or from scratch)
here's what I'm proposing for the Python interface:

B=trimult(A, B, uplo='U', side='L', inplace=True)
B=trisolve(A, B, uplo, side, inplace)

A is triangular, B is rectangular or a vector. We could trim the uplo and
side arguments if there were 'triangular matrix' classes that kept track of
whether they're upper or lower, but that would be more trouble than it's
worth, no?

I'm thinking if inplace=True, B will only be overwritten if it's possible to
recycle the memory efficiently, so the only safe usage is to use the return
value.

How does this sound?

Anand
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/scipy-dev/attachments/20071128/46e2b21a/attachment.html>