Can we use scipy to do faster LU decomposition for band matrices?
Lapack's *gbsv routine computes the LU decomp of an input banded matrix. 
From python, you can use either its f2py wrapper (see e.g. the source of scipy.linalg.solve_banded for example usage) or drop to Cython and use scipy.linalg.cython_lapack bindings.
We know that elimination requires roughly 1/3 n^3 operations, and if we use LU decomposition stored in memory, it is reduced to n^2 operations. If we have a band matrix with w upper and lower diagonals, we can skip the zeros and bring it down to about nw^2 operations, and if we use LU decomposition, it can be done in about 2nw operations.
In scipy.linalg, we have lu_factor and lu_solve, but they do not seem to be optimized for band matrices. We also have solve_banded, but it directly solves Ax=b. How can we do an efficient LU decomposition for banded matrices and efficiently perform forward and backward elimination with banded triangular L and U?