fix promotion in bidiagonal backslash, At_ldiv_B, and Ac_ldiv_B. add missing bidiagonal matrix conversions.

[Sacha0]

This snippet

m = 10
intvec = fill(1, m);
lbdintmat = Bidiagonal(fill(2, m), fill(1, m-1), false);
(\)(lbdintmat, intvec)

yields

ERROR: InexactError()
 in naivesub! at linalg/bidiag.jl:253
 [inlined code] from linalg/bidiag.jl:206
 in \ at linalg/bidiag.jl:267
 in eval at ./boot.jl:265

due to promotion in the relevant backslash definition not accounting for usually-necessary divisions in forward/back substitution. This pull request fixes the promotion.

Needs tests. Question on that front. The linalg/bidiag tests don't catch this by nature of testing only floating-point types. Two potential approaches:

1. Minimal modification. Add targeted tests for this issue.
2. Significant modification. Expand the set and combinations of types tested.

? Thanks, and best!

[Sacha0]

Seems Appveyor timed out on x86_64 here as well?

[Sacha0]

Expanding the linalg/bidiag tests to integer types revealed the same issue with A(t|c)_ldiv_B for integer bidiagonal A and integer matrix B: The last two calls in this snippet

m = 10
intmat = fill(1, m, m);
lbdintmat = Bidiagonal(fill(2, m), fill(1, m-1), false);
At_ldiv_B(lbdintmat, intmat)
Ac_ldiv_B(lbdintmat, intmat)

both fail as above. The offending lines are

julia/base/linalg/bidiag.jl

Lines 240 to 241 in 25c3659

	Ac_ldiv_B(A::Union{Bidiagonal, AbstractTriangular}, B::AbstractMatrix) = Ac_ldiv_B!(A,copy(B))
	At_ldiv_B(A::Union{Bidiagonal, AbstractTriangular}, B::AbstractMatrix) = At_ldiv_B!(A,copy(B))

. The same issue exists with a triangular matrix in place of the bidiagonal, but not where B is a StridedVecOrMat due to shadowing by more specialized methods in base/linalg/triangular.jl. Sparse B demonstrates the issue though:

m = 10
intmat = fill(1, m, m);
ltintmat = LowerTriangular(rand(1:10, m, m));
At_ldiv_B(ltintmat, sparse(intmat))
Ac_ldiv_B(ltintmat, sparse(intmat))

So the two lines linked above need fixing as well.

The linalg/triangular tests fail to catch this issue for triangular matrices as the inner loop over eltypB (

julia/test/linalg/triangular.jl

Line 299 in 25c3659

for eltyB in (Float32, Float64, Complex64, Complex128)

) contains only real and complex floating-point types, and they test only StridedVecOrMat Bs in any case. That statement holds for the expanded tests in #14508 as well. So on top of #14508, perhaps the eltypB loop in test/linalg/triangular.jl should expand to some integer types, and a test specific to this issue should go into the sparse tests.

I will fix the offending lines and update this PR.

[Sacha0]

Should be ready for review. A few points:

The promotion fixes and expanded tests require the added convertmethods.

The promotion issue involves a generic linear algebra subtlety, described below. Due to that subtlety, the fix for the triangular case is more extensive and largely orthogonal to this PR. So instead of pursuing that fix here, I decoupled the two cases and added a fix-me note to the triangular code. If the approach to the fix in this PR / described below receives approval, I will open a separate issue/PR for the triangular case.

The generic linear algebra subtlety, closely related to recent discussion in JuliaLang/LinearAlgebra.jl#294: To fix the promotion issue's bidiagonal case, I replaced

Ac_ldiv_B(A::Union{Bidiagonal, AbstractTriangular}, B::AbstractMatrix) = Ac_ldiv_B!(A,copy(B))
At_ldiv_B(A::Union{Bidiagonal, AbstractTriangular}, B::AbstractMatrix) = At_ldiv_B!(A,copy(B))

and the related \ method, with

for (f,g) in ((:\, :A_ldiv_B!), (:At_ldiv_B, :At_ldiv_B!), (:Ac_ldiv_B, :Ac_ldiv_B!))
    @eval begin
        function ($f){TA<:Number,TB<:Number}(A::Bidiagonal{TA}, B::AbstractVecOrMat{TB})
            TAB = typeof((zero(TA)*zero(TB) + zero(TA)*zero(TB))/one(TA))
            ($g)(convert(AbstractArray{TAB}, A), copy_oftype(B, TAB))
        end
        ($f)(A::Bidiagonal, B::AbstractVecOrMat) = ($g)(A, copy(B))
    end
end

These definitions differ slightly from the related generic promotion method definitions in base/linalg/triangular.jl: I qualified TA and TB in the promotion methods with <:Number and added more general fallback methods. Reason: The generic promotion methods in base/linalg/triangular.jl will fail for matrices with element type T lacking zero(T)/one(T) methods, for example block-matrix constructions; this problem necessitates a more general fallback method. The <:Number qualification prevents this more general fallback method from colliding with the generic promotion method; <:Number seems the most general qualification that makes sense(?), though this is a great use-case for traits (dispatch on hasone(T), haszero(T) or so).

If this seems cogent, I will open an issue/PR to carry these changes to base/linalg/triangular.jl. Thanks, and best!

[tkelman]

is this backwards?

[Sacha0]

Which part? Those lines follow

julia/base/linalg/triangular.jl

Lines 23 to 24 in 9bab4da

	convert{Tnew,Told,S}(::Type{$t{Tnew}}, A::$t{Told,S}) = (Anew = convert(AbstractMatrix{Tnew}, A.data); $t(Anew))
	convert{Tnew,Told,S}(::Type{AbstractMatrix{Tnew}}, A::$t{Told,S}) = convert($t{Tnew}, A)

which seem the same? Also, tests pass which are dependent on the functionality?

[Sacha0]

Oh, do you mean the comment? Is the AbstractMatrix{Tnew}(->Bidiagonal{Tnew}) part unclear?

[tkelman]

maybe separate the 2 comments and use words for what the -> part is indicating

[Sacha0]

Shall I just nix the comment altogether? Those lines seem more or less self-explanatory. Or do you prefer the comment? If the latter, I'll break it in two. Thanks!

[tkelman]

the first line is mostly clear, converting the data fields, but the second line isn't quite so obvious. converting to an abstract type is a little bit ambiguous

[andreasnoack]

We use these kinds of conversions for element type promotion for all the special matrices. Explaining this in a comment would probably be useful.

[Sacha0]

Revised accordingly. Thanks!

[Sacha0]

Is this in shape to merge? If so I will open an issue to track the triangular case. Fixing the triangular case appears more involved than fixing the bidiagonal case, and might be better to address after JuliaLang/LinearAlgebra.jl#283 (which reminds me that I need clean up #14447.) Thanks, and best!

[andreasnoack]

Did you check if the methods we discuss were shadowed?

[Sacha0]

Did you check if the methods we discuss were shadowed?

The general promotion methods in linalg/triangular.jl only shadow [the two fallbacks moved from linalg/bidiag.jl] for StridedVecOrMat right-hand-sides. To illustrate, the fallbacks handle

ltA = LowerTriangular(eye(5) + randn(5, 5)/10)
svb = sparsevec([1], [1.0], 5)
At_ldiv_B(ltA, svb)

Hence I believe the least dangerous / potentially disruptive approach is to leave those fallbacks in place till the triangular case is fixed properly. Thoughts? Thanks, and best!

[andreasnoack]

What if eltype(B) is not a Number and not able to store the solution? On the top of my head, that could be something like eltype(A)=Matrix{Float64} and eltype(B)=Matrix{Int}

[Sacha0]

Then the fallback fails, as would the promotion method were the Number restriction relaxed. Graceful handling of block matrices seems a broader issue. I've been thinking about that issue in the background following #14490 and hope to partially address it, but I am not prepared to begin doing so here; rather, at some point I will begin with a followup to #14490 and #14493. But perhaps I miss a simple and graceful patch? Thanks again!

[andreasnoack]

Okay. That's fine. This still covers move cases than it used to.

[Sacha0]

Thanks for the thoughtful review and merge!