[SciPy-dev] Scipy-dev Digest, Vol 49, Issue 14

[Nathan Bell]

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.

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.

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

[1] http://en.wikipedia.org/wiki/Gauss-Seidel_method
[2] http://projects.scipy.org/scipy/scipy/browser/trunk/scipy/sandbox/multigrid/multigridtools/relaxation.h

> Yeah, if I'm going to code this triangular stuff anyway I may as well do it
> in such a way that's useful to other people.
>
> Any thoughts on wrapping a library vs. writing from scratch?

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.  The current
bindings automatically upcast their arguments to the appropriate type,
which simplifies the Python side a bit.  Special functions to sort the
column/row indices of a CSR/CSC matrix may not exist in other
libraries.  Other libraries may not support numpy's integer of complex
types out of the box.  I didn't trust that other libraries had
efficient implementations of level 3 operations between CSR/CSC types,
conversions among the formats, etc.

Ultimately, I think making the backend code amenable to SWIG and
scipy.sparse saved more effort than starting with an existing library.

-- 
Nathan Bell wnbell at gmail.com