Solve sparse upper triangular system
Question: I'm trying to figure out how to efficiently solve a sparse triangular system,
Au*x = b in scipy sparse.

For example, we can construct a sparse upper triangular matrix, Au, and a
right hand side b with:

    
    
    import scipy.sparse as sp
    import scipy.sparse.linalg as sla
    import numpy as np
    
    n = 2000
    A = sp.rand(n, n, density=0.4) + sp.eye(n)
    Au = sp.triu(A).tocsr()
    b = np.random.normal(size=(n))
    

We can get a solution to the problem using spsolve, however it is clear that
the triangular structure is not being taken advantage of. This can be
demonstrated by timing the solution and comparing it to the solve method in
splu. (Here using iPython's %time magic)

    
    
    %time x1 = sla.spsolve(Au,b)
    CPU times: user 3.63 s, sys: 79.1 ms, total: 3.71 s
    Wall time: 1.1 s
    
    %time Au_lu = sla.splu(Au)
    CPU times: user 3.61 s, sys: 62.2 ms, total: 3.67 s
    Wall time: 1.08 s
    
    %time x2 = Au_lu.solve(b)
    CPU times: user 25 ms, sys: 332 µs, total: 25.4 ms
    Wall time: 7.01 ms
    

As Au is already upper-triangular the call to splu really shouldn't be doing
much of anything, however, as n gets large this call becomes prohibitively
expensive (as does the use of spsolve), while the solve time remains small.

Is there any way of using superLU's triangular solver without first calling
splu? Or is there a better way of doing this altogether?


Answer: I'm afraid this isn't terribly instructive, but have you tried changing the
column permutation? When I use "NATURAL", I get huge speedups.

    
    
    %time x1 = sla.spsolve(Au, b, permc_spec="NATURAL")
    CPU time: user 46.7 ms, sys: 0 ns, total: 46.7 ms
    Wall time: 49 ms
    

For me, it's not quite as fast as using the solve methods of the splu function
output, but it seems to get reasonably close (and avoids calling splu).
Perhaps that will suffice? [Scipy
Docs](http://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.linalg.spsolve.html#scipy.sparse.linalg.spsolve)

