JuliaLang/julia — pull request #1243 by timholy (2012-08-31 14:09:25+00:00)
Title: RFC: Improved tridiagonal and new Woodbury matrix algebra

This commit:
- Converts Tridiagonal to being Lapack-compatible
- Wraps Lapack's major tridiagonal routines
- Adds a new "Woodbury" type for solving equations using the Woodbury matrix identity
- Provides a set of tests for this functionality

Issues that might be worth discussing before merging this:
- Should the tridiagonal routines be moved into ``linalg_lapack.jl`` and friends? Conceptually, tridiagonal routines are sparse, which is why I've been putting them in ``extras/sparse.jl``. However,
* Inside ``linalg_lapack.jl``, at least some "banded" matrix routines, like ``gbtrf``, are Julia-wrapped. Tridiagonal matrices are simply an extreme case of banded matrices. (As far as I can tell, there's no way to actually _use_ these banded matrix routines currently, because the basic Lapack wrappers are not exported and there is no higher-level routine that calls them.)
* There's some awkwardness to re-loading the OpenBLAS library, particular if the user has chosen a different BLAS library
- Currently I don't provide an ``lu`` wrapper for LU-decomposition, but there is a "private" (underscore) wrapper for tridiagonal LU decomposition, a wrapper to use the LU decomposition to solve equations, and a test of this functionality in ``test/sparse.jl``. The reason I don't provide a nice ``lu`` function is that the LU-decomposition of a tridiagonal matrix requires the second off-diagonal. This can actually be stored in the Tridiagonal type (and I do so), but the meaning of the object is not what it once was. There are several possible ways around this:
* Have ``lu(M::Tridiagonal)`` return full (i.e., dense) ``L`` and ``U`` matrices. This stinks so bad, it is hardly worth even mentioning.
* Have ``lu(M::Tridiagonal)`` do what the underscore Lapack wrapper does now: return an object of type Tridiagonal, but this is something of a lie, because the second off-diagonal stored in the ``dutmp`` field is actually meaningful (it just happens to be convenient storage that is already allocated). The second output would be ``ipiv`` or a permutation constructed from it. In other words, the return-value syntax of ``lu`` would differ from that of dense LU-decomposition.
* Create BandedMatrix types in Julia. This may be the most ambitious route, with the advantage that the current return syntax of ``lu`` could be preserved. However, there is a performance penalty here, because the LU decomposition of a tridiagonal can be stored and manipulated more efficiently than what is required for general banded matrix support.
* Introduce a new alternative syntax ``lu(LU::LUMatrixDecomposition, A::Matrix)``, a function ``solve(LU::LUMatrixDecomposition, rhs::StridedVecOrMat)``, and types to support this. This syntax allows you to control the output format of the result, and uses pre-allocated storage consistent with issue #1115. It would also allow for maximally-efficient usage of the Lapack routines, without the need to call ``tril`` and ``triu`` on the result. That said, compared to the ``O(N^3)`` cost of LU-decomposition, the ``O(N^2)`` overhead of ``tril`` and ``triu`` should be pretty inconsequential. So the main advantage of this route is simply that you can customize the return type appropriately to the input type, which is something you can't do with the current ``lu`` syntax.

For this second issue, I personally think the fourth option is the best. Interested in hearing feedback from others.