[SciPy-User] iterative matrix methods seem slow

[Richard Sharp]

Hello,

I'm attempting to solve a largish sparse, banded diagonal, linear
system resulting from a finite difference approximation of a
diffusion-like equation.

I can use `spsolve` on smaller matrices, but get an odd memory error
when they get to around (640k)^2:

...Can't expand MemType 1

I've also looked into using some of the iterative solvers, but they
are all painfully slow.  For example, on a (40k)^2 system `spsolve`
runs in 0.8s while gmres takes 4.5s  Bumping up to (90k)^2 takes 2.4s
with `spsolve` and 15.4s on `gmres`.  The other methods don't work, or
run so long I don't wait for them to converge.

I've tried using a Jacobi preconditioner and making good
approximations to the solution for `gmres`, but I only end up with
slower times.  I think I'm doing something wrong, because It's my
impression that the iterative methods should run pretty quickly
especially compared against a direct solver.

My code looks something like this:

    #A holds diagonals
    matrix = spdiags(A,[-m,-1,0,1,m],n*m,n*m,"csr")
    if direct:
        result = spsolve(matrix,b)
    else:
        print ' system too large trying iteration'
        result = scipy.sparse.linalg.gmres(matrix,b)[0]

I'd appreciate any help with anything I could be doing wrong with
setting up the system, making the calls, or a fundamental
misunderstanding of the methods.

Thanks for any help,
Rich

--
Richard P. Sharp Jr.
Lead Software Developer
Natural Capital Project
Stanford University, U Minnesota, TNC, WWF
371 Serra Mall
Stanford, CA 94305
http://www.stanford.edu/~rpsharp/