[SciPy-User] solving large sparse linear system with Laplacian matrix
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
> 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
> 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:
Convert a symmetric banded matrix to symmetric tridiagonal form:
The tridiagonal form is useful because you can then use this
specialized routine for solving symmetric tridiagonal systems:
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
Hope this helps,
More information about the SciPy-User