Merge pull request #14506 from Sacha0/bidiagbackslashprom

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

Author: andreasnoack

base/linalg/bidiag.jl

@@ -77,6 +77,13 @@ function convert{T}(::Type{Tridiagonal{T}}, A::Bidiagonal)
 end
 promote_rule{T,S}(::Type{Tridiagonal{T}}, ::Type{Bidiagonal{S}})=Tridiagonal{promote_type(T,S)}
 
+# No-op for trivial conversion Bidiagonal{T} -> Bidiagonal{T}
+convert{T}(::Type{Bidiagonal{T}}, A::Bidiagonal{T}) = A
+# Convert Bidiagonal{Told} to Bidiagonal{Tnew} by constructing a new instance with converted elements
+convert{Tnew,Told}(::Type{Bidiagonal{Tnew}}, A::Bidiagonal{Told}) = Bidiagonal(convert(Vector{Tnew}, A.dv), convert(Vector{Tnew}, A.ev), A.isupper)
+# When asked to convert Bidiagonal{Told} to AbstractMatrix{Tnew}, preserve structure by converting to Bidiagonal{Tnew} <: AbstractMatrix{Tnew}
+convert{Tnew,Told}(::Type{AbstractMatrix{Tnew}}, A::Bidiagonal{Told}) = convert(Bidiagonal{Tnew}, A)
+
 big(B::Bidiagonal) = Bidiagonal(big(B.dv), big(B.ev), B.isupper)
 
 ###################
@@ -230,7 +237,6 @@ function A_ldiv_B!(A::Union{Bidiagonal, AbstractTriangular}, B::AbstractMatrix)
     end
     B
 end
-
 for func in (:Ac_ldiv_B!, :At_ldiv_B!)
     @eval function ($func)(A::Union{Bidiagonal, AbstractTriangular}, B::AbstractMatrix)
         nA,mA = size(A)
@@ -247,9 +253,6 @@ for func in (:Ac_ldiv_B!, :At_ldiv_B!)
         B
     end
 end
-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))
-
 #Generic solver using naive substitution
 function naivesub!{T}(A::Bidiagonal{T}, b::AbstractVector, x::AbstractVector = b)
     N = size(A, 2)
@@ -272,9 +275,15 @@ function naivesub!{T}(A::Bidiagonal{T}, b::AbstractVector, x::AbstractVector = b
     x
 end
 
-function \{T,S}(A::Bidiagonal{T}, B::AbstractVecOrMat{S})
-    TS = typeof(zero(T)*zero(S) + zero(T)*zero(S))
-    TS == S ? A_ldiv_B!(A, copy(B)) : A_ldiv_B!(A, convert(AbstractArray{TS}, B))
+### Generic promotion methods and fallbacks
+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
 
 factorize(A::Bidiagonal) = A

base/linalg/triangular.jl

@@ -1366,6 +1366,11 @@ for (f, g) in ((:/, :A_rdiv_B!), (:A_rdiv_Bc, :A_rdiv_Bc!), (:A_rdiv_Bt, :A_rdiv
         end
     end
 end
+### Fallbacks brought in from linalg/bidiag.jl while fixing #14506.
+# Eventually the above promotion methods should be generalized as
+# was done for bidiagonal matrices in #14506.
+At_ldiv_B(A::AbstractTriangular, B::AbstractVecOrMat) = At_ldiv_B!(A, copy(B))
+Ac_ldiv_B(A::AbstractTriangular, B::AbstractVecOrMat) = Ac_ldiv_B!(A, copy(B))
 
 # Complex matrix logarithm for the upper triangular factor, see:
 #   Al-Mohy and Higham, "Improved inverse  scaling and squaring algorithms for

test/linalg/bidiag.jl

@@ -8,16 +8,28 @@ n = 10 #Size of test matrix
 srand(1)
 
 debug && println("Bidiagonal matrices")
-for relty in (Float32, Float64, BigFloat), elty in (relty, Complex{relty})
+for relty in (Int, Float32, Float64, BigFloat), elty in (relty, Complex{relty})
     debug && println("elty is $(elty), relty is $(relty)")
-    dv = convert(Vector{elty}, randn(n))
-    ev = convert(Vector{elty}, randn(n-1))
-    b = convert(Matrix{elty}, randn(n, 2))
-    c = convert(Matrix{elty}, randn(n, n))
-    if (elty <: Complex)
-        dv += im*convert(Vector{elty}, randn(n))
-        ev += im*convert(Vector{elty}, randn(n-1))
-        b += im*convert(Matrix{elty}, randn(n, 2))
+    if relty <: AbstractFloat
+        dv = convert(Vector{elty}, randn(n))
+        ev = convert(Vector{elty}, randn(n-1))
+        b = convert(Matrix{elty}, randn(n, 2))
+        c = convert(Matrix{elty}, randn(n, n))
+        if (elty <: Complex)
+            dv += im*convert(Vector{elty}, randn(n))
+            ev += im*convert(Vector{elty}, randn(n-1))
+            b += im*convert(Matrix{elty}, randn(n, 2))
+        end
+    elseif relty <: Integer
+        dv = convert(Vector{elty}, rand(1:10, n))
+        ev = convert(Vector{elty}, rand(1:10, n-1))
+        b = convert(Matrix{elty}, rand(1:10, n, 2))
+        c = convert(Matrix{elty}, rand(1:10, n, n))
+        if (elty <: Complex)
+            dv += im*convert(Vector{elty}, rand(1:10, n))
+            ev += im*convert(Vector{elty}, rand(1:10, n-1))
+            b += im*convert(Matrix{elty}, rand(1:10, n, 2))
+        end
     end
 
     debug && println("Test constructors")
@@ -89,23 +101,24 @@ for relty in (Float32, Float64, BigFloat), elty in (relty, Complex{relty})
         condT = cond(map(Complex128,Tfull))
         x = T \ b
         tx = Tfull \ b
+        promty = typeof((zero(relty)*zero(relty) + zero(relty)*zero(relty))/one(relty))
         @test_throws DimensionMismatch Base.LinAlg.naivesub!(T,ones(elty,n+1))
-        @test norm(x-tx,Inf) <= 4*condT*max(eps()*norm(tx,Inf), eps(relty)*norm(x,Inf))
+        @test norm(x-tx,Inf) <= 4*condT*max(eps()*norm(tx,Inf), eps(promty)*norm(x,Inf))
         @test_throws DimensionMismatch T \ ones(elty,n+1,2)
         @test_throws DimensionMismatch T.' \ ones(elty,n+1,2)
         @test_throws DimensionMismatch T' \ ones(elty,n+1,2)
         if relty != BigFloat
             x = T.'\c.'
             tx = Tfull.' \ c.'
-            @test norm(x-tx,Inf) <= 4*condT*max(eps()*norm(tx,Inf), eps(relty)*norm(x,Inf))
+            elty <: AbstractFloat && @test norm(x-tx,Inf) <= 4*condT*max(eps()*norm(tx,Inf), eps(promty)*norm(x,Inf))
             @test_throws DimensionMismatch T.'\b.'
             x = T'\c.'
             tx = Tfull' \ c.'
-            @test norm(x-tx,Inf) <= 4*condT*max(eps()*norm(tx,Inf), eps(relty)*norm(x,Inf))
+            @test norm(x-tx,Inf) <= 4*condT*max(eps()*norm(tx,Inf), eps(promty)*norm(x,Inf))
             @test_throws DimensionMismatch T'\b.'
             x = T\c.'
             tx = Tfull\c.'
-            @test norm(x-tx,Inf) <= 4*condT*max(eps()*norm(tx,Inf), eps(relty)*norm(x,Inf))
+            @test norm(x-tx,Inf) <= 4*condT*max(eps()*norm(tx,Inf), eps(promty)*norm(x,Inf))
             @test_throws DimensionMismatch T\b.'
         end
 
@@ -133,11 +146,13 @@ for relty in (Float32, Float64, BigFloat), elty in (relty, Complex{relty})
         @test_throws ArgumentError diag(T,n+1)
 
         debug && println("Eigensystems")
-        d1, v1 = eig(T)
-        d2, v2 = eig(map(elty<:Complex ? Complex128 : Float64,Tfull))
-        @test_approx_eq isupper?d1:reverse(d1) d2
-        if elty <: Real
-            Test.test_approx_eq_modphase(v1, isupper?v2:v2[:,n:-1:1])
+        if relty <: AbstractFloat
+            d1, v1 = eig(T)
+            d2, v2 = eig(map(elty<:Complex ? Complex128 : Float64,Tfull))
+            @test_approx_eq isupper?d1:reverse(d1) d2
+            if elty <: Real
+                Test.test_approx_eq_modphase(v1, isupper?v2:v2[:,n:-1:1])
+            end
         end
 
         debug && println("Singular systems")
@@ -161,8 +176,8 @@ for relty in (Float32, Float64, BigFloat), elty in (relty, Complex{relty})
         @test convert(elty,-1.0) * T == Bidiagonal(-T.dv,-T.ev,T.isupper)
         @test T * convert(elty,-1.0) == Bidiagonal(-T.dv,-T.ev,T.isupper)
         for isupper2 in (true, false)
-            dv = convert(Vector{elty}, randn(n))
-            ev = convert(Vector{elty}, randn(n-1))
+            dv = convert(Vector{elty}, relty <: AbstractFloat ? randn(n) : rand(1:10, n))
+            ev = convert(Vector{elty}, relty <: AbstractFloat ? randn(n-1) : rand(1:10, n-1))
             T2 = Bidiagonal(dv, ev, isupper2)
             Tfull2 = full(T2)
             for op in (+, -, *)