# [SciPy-User] solving large sparse linear system with Laplacian matrix

David Warde-Farley dwf@cs.toronto....
Fri Oct 30 04:28:13 CDT 2009

```On 30-Oct-09, at 3:46 AM, Emmanuelle Gouillart wrote:

> I'm looking for a memory-savvy algorithm in scipy.sparse.linalg,
> that I
> could use to solve a very large sparse linear system. The system can
> be
> written as AX = B, where A is a symmetric band matrix with non-zero
> elements on only 4 or 6 upper (and lower) diagonals. The shape of A is
> NxN where N is very large (up to 10e6, if possible...), the shape of B
> and X is NxM, where M is much smaller (up to 50, say).
> Unfortunately, A
> is symmetric but not positive-definite as the sum of each row and
> column
> is null (A is a Laplacian matrix). B is also sparse (it is actually a
> block of the original Laplacian matrix).

I'm not sure they exist in scipy.sparse.linalg, but looking into your
problem a little further it appears that several options exist in
LAPACK (but are not wrapped in SciPy, yet). Conveniently, they have
built-in support for multiple right-hand-sides.

Solve a general banded system of linear equations:
http://www.netlib.org/lapack/double/dgbsv.f

Convert a symmetric banded matrix to symmetric tridiagonal form:
http://www.netlib.org/lapack/double/dsbtrd.f

The tridiagonal form is useful because you can then use this
specialized routine for solving symmetric tridiagonal systems:
http://www.netlib.org/lapack/double/dptsv.f

So, if both memory and performance are at issue, I'd recommend looking
into wrapping dsbtrd and dptsv with f2py (these are for doubles, hence
the 'd' - if you're interested in single precision use the 's'
versions). Note that they require you pass in the band matrix in a
particular compressed format, which is kid of a pain, but on the other
hand _will_ save you space over CSR or anything like that. The layout
you need to use for symmetric band matrices is near the bottom of this
page:

http://www.netlib.org/lapack/lug/node124.html

Hope this helps,

David
```