allow A_ldiv_B! specialization in generic triangular solves

[Sacha0]

In 53978e7/base/linalg/bidiag.jl#L206-L208, generic methods for A(t|c)_ldiv_B! directly call naivesub! for A_ldiv_B! functionality:

A_ldiv_B!(A::Union{Bidiagonal, AbstractTriangular}, b::AbstractVector) = naivesub!(A, b)
At_ldiv_B!(A::Union{Bidiagonal, AbstractTriangular}, b::AbstractVector) = naivesub!(transpose(A), b)
Ac_ldiv_B!(A::Union{Bidiagonal, AbstractTriangular}, b::AbstractVector) = naivesub!(ctranspose(A), b)

This PR simply replaces the relevant naivesub! calls with A_ldiv_B!,

A_ldiv_B!(A::Union{Bidiagonal, AbstractTriangular}, b::AbstractVector) = naivesub!(A, b)
At_ldiv_B!(A::Union{Bidiagonal, AbstractTriangular}, b::AbstractVector) = A_ldiv_B!(transpose(A), b)
Ac_ldiv_B!(A::Union{Bidiagonal, AbstractTriangular}, b::AbstractVector) = A_ldiv_B!(ctranspose(A), b)

which gives dispatch a chance to find better methods than naivesub!. Impact:

# setup
ml, ms = 10000, 1000;
bo, Bo = rand(ml), rand(ms, ms);
As = eye(ms) + randn(ms, ms)/(2ms);
Al = eye(ml) + randn(ml, ml)/(2ml);
ltAs = LowerTriangular(As);
ltAl = LowerTriangular(Al);

Without PR:

# warmup hidden
julia> B = copy(Bo);
julia> @time Base.LinAlg.At_ldiv_B!(ltAs, B);
    8.259813 seconds (6.00 k allocations: 7.451 GB, 10.34% gc time)
julia> b = copy(bo);
julia> @time Base.LinAlg.At_ldiv_B!(ltAl, b);
    1.412599 seconds (10 allocations: 762.940 MB, 4.28% gc time)

With PR:

# warmup hidden
julia> B = copy(Bo);
julia> @time Base.LinAlg.At_ldiv_B!(ltAs, B);
    5.691988 seconds (7.00 k allocations: 7.451 GB, 14.62% gc time)
julia> b = copy(bo);
julia> @time Base.LinAlg.At_ldiv_B!(ltAl, b);
    0.691178 seconds (57 allocations: 762.942 MB, 0.53% gc time)

Best!

[yuyichao]

AppVeyor Win64 failure is #14395 (I aborted it after no update in output for 5mins). At least this time the failing malloc is among the ones we check and the error message is quite clear what's happening.

[andreasnoack]

Good catch, but I'm wondering if we should define the missing methods instead of computing the transpose. Ac_ldiv_B! is defined for BlasFloats, but not At_ldiv_B!. This should be easy to add, but we lack the generic fallbacks for both versions.

[Sacha0]

I'm wondering if we should define the missing methods instead of computing the transpose. Ac_ldiv_B! is defined for BlasFloats, but not At_ldiv_B!. This should be easy to add, but we lack the generic fallbacks for both versions.

My thinking as well. If you would welcome more significant changes (like these), I would be happy to prepare a few related pull requests (beyond what we discussed in JuliaLang/LinearAlgebra.jl#283.) Being new I've hesitated in submitting more sweeping modifications uninvited / undiscussed.

[andreasnoack]

It would be great if you could prepare pull requests with the missing methods.

[andreasnoack]

This is an improvement over what we have. So I'll merge, but specialized methods for transposed bidiagonal and triangular solves would be good to have.