[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


More information about the SciPy-User mailing list