sympy/sympy — pull request #8083 by pbrady (2014-09-23 15:36:37+00:00)
Title: Add banded diagonal storage and solver

Finiite-difference and finite-volume schemes tend to yield matrices with a banded structure. For such matrices, the DOK implementation is far from optimal so this PR implements a matrix diagonal storage form and solver.

Sample session:

```python
$ ipython
Python 3.4.1 |Continuum Analytics, Inc.| (default, Sep 2 2014, 14:00:37)
Type "copyright", "credits" or "license" for more information.

IPython 2.2.0 -- An enhanced Interactive Python.
Anaconda is brought to you by Continuum Analytics.
Please check out: http://continuum.io/thanks and https://binstar.org
? -> Introduction and overview of IPython's features.
%quickref -> Quick reference.
help -> Python's own help system.
object? -> Details about 'object', use 'object??' for extra details.

In [1]: from sympy import SparseMatrix, Matrix, randMatrix

In [2]: from numpy.random import randint

In [3]: from sympy.matrices import sparsetools

In [4]: np = 100

In [5]: sol = randMatrix(np, 5)

In [6]: upper, lower = 1, 3 # bandwidths

In [7]: def fd_mat(i, j, upper=upper, lower=lower, np=np):
for shift in range(-upper, lower+1):
if i==j+shift:
return randint(shift, high=2*np//(abs(shift)+1))
...:

In [8]: s = SparseMatrix(np, np, fd_mat)

In [9]: (sparsetools.banded_LUsolve(s, sol) - s.LUsolve(sol)).is_zero
Out[9]: True

In [10]: %timeit sparsetools.banded_LUsolve(s, sol)
1 loops, best of 3: 976 ms per loop

In [11]: %timeit s.LUsolve(sol)
1 loops, best of 3: 14.6 s per loop
```

The new routine `banded_LUsolve` can simply take a DOK sparse matrix, convert it into diag storage, compute the LU and solve. However, if one is going to be solving the same system many times with different rhs it can be beneficial to precompute the LU. This is available as well (somewhat hidden because it's an 'expert' option)

```python
In [12]: l_and_u, LU = sparsetools._banded_LUdecomposition(s)

In [13]: (sparsetools.banded_LUsolve(LU, sol, l_and_u, have_LU=True) - s.LUsolve(sol)).is_zero
Out[13]: True

In [14]: %timeit sparsetools.banded_LUsolve(LU, sol, l_and_u, have_LU=True)
1 loops, best of 3: 829 ms per loop
```

Here are items that need to be addressed to move forward

* [ ] Doctests?
* [ ] Is sparse_tools the appropriate place?
* [ ] tests

@tjl do you want to test this out?