Add wrappers for LAPACK pbtrf, pbtrs

[nignatiadis]

This wraps the LAPACK routines for the Cholesky factorization of symmetric (or Hermitian) positive definite banded matrices.

[tkelman]

This would probably be better for https://github.com/ApproxFun/BandedMatrices.jl as long as we don't have types for banded matrices in base.

[nignatiadis]

Thanks for the quick reply and the link! I somehow missed this package.

I submitted this pull request, mainly because the functionality was very similar to gbtrf! and gbtrs! (LU decomposition for banded matrices), which are in Base.

[tkelman]

The banded LU wrappers have been in base for some time, before BandedMatrices.jl was registered. It may make sense to move them out too.

[nignatiadis]

Yes, makes sense. I'll try opening a pull request/issue there too. It does not directly fit there right now either, since these matrices are not only banded, but also symmetric and positive definite (and currently there are no corresponding types in BandedMatrices.jl). I guess https://github.com/JuliaStats/PDMats.jl could also work.

[andreasnoack]

I think it would make sense to try to define the relevant types in https://github.com/ApproxFun/BandedMatrices.jl but the maintainers of that package will be able to answer that question better. Personally, I don't think PDMats is the right choice very often because, in computation, when you want to exploit PD then you almost always use the Cholesky. Therefore having a type for positive definiteness seems redundant to me.

[gasagna]

In BandedMatrices.jl we only use the gb* functions for lu factorisation of generic banded matrices. This code would be welcome. However, currently we do not have specialised types but it would not be a problem to extend the type hierarchy to pd matrices. @dlfivefifty

edit: in that package we call the Julia wrappers to gb* which are already in LinAlg. I would not mind if we could do the same for pb* from BandedMatrices.jl. Having all lapack functions wrapped in base is convenient for package developers as it avoid having the same implementation scattered in too many places.

[dlfivefifty]

I'd be happy to have

PDBandedMatrix{T} <: AbstractBandedMatrix{T}

in BandedMatrices.jl.

It might also make sense to move it to a new organization, say, JuliaMatrices? This could also incorporate the ToeplitzMatrices.jl package.

PS Anyone here have opinions on "BandedMatrix" vs "BandMatrix"? I'm considering renaming it because "BandMatrix" is shorter and is also a common term, but I can't decide.

[tkelman]

I have a slight preference for Banded as it's an adjective, Band is a noun.

Might be worth trying to use the base Symmetric and Hermitian wrappers, since symmetry is the important part from a data structure perspective. Positive definiteness is important for behavior and how many algorithms will behave, but it's not always verifiable by construction until you try a Cholesky factorization on it. There aren't as many routines in lapack for symmetric indefinite banded matrices, but they come up in a fair number of applications.

Having all lapack functions wrapped in base is convenient for package developers as it avoid having the same implementation scattered in too many places.

I can see where you're coming from, but we only wrap a subset of lapack, gmp, mpfr etc, which are all huge libraries and the base wrappers don't pretend to be comprehensive. If code is unused in base and has no technical reason to be here, it's better suited for a package IMO.

[gasagna]

@nignatiadis I think a good start would be have julia's gbrtf/s code and the p* versions in a lapack.jl file in BandedMatrices.jl. We can then have a cholfact(::BandedMatrix) method, similarly to the lufact method that exists there (largely inspired by julia base code).

Currently, for lufact we copy the bands into a larger storage with 2l+u+1 rows (l/u -> lower/upper bandwidths) required to accommodate the fill-in due to factorisation. You could do the same. Then, when we reorganise how the internal storage is handled we can add custom types, and exploit the structure more efficiently.

@tkelman Has anyone ever considered Lapack.jl package?

[nignatiadis]

Thanks to everyone for the helpful and instructive replies! I agree that it would make a lot of sense to have a LAPACK.jl package, which wraps all LAPACK functions not in Base. Else I suspect there will be a lot of wasted, non-coordinated effort.

@gasagna I won't have any time until the weekend to work on this. The only thing I could do quickly is to just copy the content of this pull request to the blas.jl file of BandedMatrices. But maybe instead of doing that, it would make more sense to start the LAPACK.jl package and call it from BandedMatrices? We could start with gbtrf/s, pbtrf/s.

[gasagna]

@nignatiadis No worries, take your time.

I believe a LAPACK.jl package will require substantial efforts from the community. Probably it should be started and promoted by someone of the core developers.

[dlfivefifty]

👍 to LAPACK.jl This will ensure that there is a consistent order of arguments, etc.

[dlfivefifty]

This could also be moved to LAPACK.jl:

https://github.com/ApproxFun/ApproxFun.jl/blob/development/src/LinearAlgebra/hesseneigs.jl

It wraps the Hessenberg eigenvalue commands.

[kshyatt]

So do we want to add this into Base OR are people willing to make this LAPACK.jl package/have they done so?

[dlfivefifty]

I think it would be best as

JuliaLang/LAPack.jl

Or in a new organisation

JuliaMatrices/LAPack.jl

Sent from my iPhone

On 4 Jul 2016, at 09:25, Katharine Hyatt notifications@github.com wrote:

So do we want to add this into Base OR are people willing to make this LAPACK.jl package/have they done so?

—
You are receiving this because you were mentioned.
Reply to this email directly, view it on GitHub, or mute the thread.

[kshyatt]

We have a lot of things in the JuliaLang org and have been moving things out. I think a linalg focused org would be a great idea.

[dlfivefifty]

Ok would it be better for someone closer to JuliaLang to create JuliaMatrices? Otherwise I can do it and add whomever wants admin access.

(It's probably better if I'm not the primary admin.)

Sent from my iPhone

On 4 Jul 2016, at 09:54, Katharine Hyatt notifications@github.com wrote:

We have a lot of things in the JuliaLang org and have been moving things out. I think a linalg focused org would be a great idea.

—
You are receiving this because you were mentioned.
Reply to this email directly, view it on GitHub, or mute the thread.

[kshyatt]

@andreasnoack seems like the obvious candidate for JuliaMatrices admin

[dlfivefifty]

👍 and he can move ToeplitzMatrices.jl there too if he's happy to.

I'd move BandedMatrices.jl in too.

Sent from my iPhone

On 4 Jul 2016, at 09:58, Katharine Hyatt notifications@github.com wrote:

@andreasnoack seems like the obvious candidate for JuliaMatrices admin

—
You are receiving this because you were mentioned.
Reply to this email directly, view it on GitHub, or mute the thread.

[andreasnoack]

https://github.com/JuliaMatrices/LAPACK.jl

[nalimilan]

JuliaMatrices sounds a bit weird for a project name hosting LAPACK.jl. Why not something more general referring to linear algebra?

[andreasnoack]

Numerical linear algebra is so closely associated with matrices that I think the name is fine.

[kshyatt]

Since the org now exists, presumably with LAPACK.jl soon or now, I'm going to close this PR. Let me know if you want a reopen.