Provides specialized Bidiagonal solvers

NEWS.md

@@ -78,8 +78,10 @@ Library improvements
       matrices of generic types ([#5255])
     - new algorithms for linear solvers, eigensystems and singular systems of `Diagonal`
       matrices of generic types ([#5263])
+    - new algorithms for linear solvers and eigensystems of `Bidiagonal`
+      matrices of generic types ([#5277])
     - specialized methods `transpose`, `ctranspose`, `istril`, `istriu` for
-      `Triangular` ([#5255])
+      `Triangular` ([#5255]) and `Bidiagonal` ([#5277])
     - new LAPACK wrappers
       - condition number estimate `cond(A::Triangular)` ([#5255])
 
@@ -123,6 +125,7 @@ Deprecated or removed
 [#5059]: https://github.com/JuliaLang/julia/issues/5059
 [#5196]: https://github.com/JuliaLang/julia/issues/5196
 [#5275]: https://github.com/JuliaLang/julia/issues/5275
+[#5277]: https://github.com/JuliaLang/julia/issues/5277
 
 Julia v0.2.0 Release Notes
 ==========================

base/linalg/bidiag.jl

@@ -1,34 +1,21 @@
-#### Specialized matrix types ####
-
-## Bidiagonal matrices
-
-
+# Bidiagonal matrices
 type Bidiagonal{T} <: AbstractMatrix{T}
-    dv::Vector{T} # diagonal
-    ev::Vector{T} # sub/super diagonal
+    dv::AbstractVector{T} # diagonal
+    ev::AbstractVector{T} # sub/super diagonal
     isupper::Bool # is upper bidiagonal (true) or lower (false)
-    function Bidiagonal{T}(dv::Vector{T}, ev::Vector{T}, isupper::Bool)
+    function Bidiagonal{T}(dv::AbstractVector{T}, ev::AbstractVector{T}, isupper::Bool)
         length(ev)==length(dv)-1 ? new(dv, ev, isupper) : throw(DimensionMismatch(""))
     end
 end
 
-Bidiagonal{T<:BlasFloat}(dv::Vector{T}, ev::Vector{T}, isupper::Bool)=Bidiagonal{T}(copy(dv), copy(ev), isupper)
-Bidiagonal{T}(dv::Vector{T}, ev::Vector{T}) = error("Did you want an upper or lower Bidiagonal? Try again with an additional true (upper) or false (lower) argument.")
+Bidiagonal{T}(dv::AbstractVector{T}, ev::AbstractVector{T}, isupper::Bool)=Bidiagonal{T}(copy(dv), copy(ev), isupper)
+Bidiagonal{T}(dv::AbstractVector{T}, ev::AbstractVector{T}) = error("Did you want an upper or lower Bidiagonal? Try again with an additional true (upper) or false (lower) argument.")
 
 #Convert from BLAS uplo flag to boolean internal
-function Bidiagonal{T<:BlasFloat}(dv::Vector{T}, ev::Vector{T}, uplo::BlasChar)
-    if uplo=='U'
-        isupper = true
-    elseif uplo=='L'
-        isupper = false
-    else
-        error("Bidiagonal can only be upper 'U' or lower 'L' but you said '$uplo''")
-    end
-    Bidiagonal{T}(copy(dv), copy(ev), isupper)
-end
+Bidiagonal(dv::AbstractVector, ev::AbstractVector, uplo::BlasChar) = Bidiagonal(copy(dv), copy(ev), uplo=='U' ? true : uplo=='L' ? false : error("Bidiagonal can only be upper 'U' or lower 'L' but you said '$uplo''"))
 
-function Bidiagonal{Td<:Number,Te<:Number}(dv::Vector{Td}, ev::Vector{Te}, isupper::Bool)
-    T = promote(Td,Te)
+function Bidiagonal{Td,Te}(dv::AbstractVector{Td}, ev::AbstractVector{Te}, isupper::Bool)
+    T = promote_type(Td,Te)
     Bidiagonal(convert(Vector{T}, dv), convert(Vector{T}, ev), isupper)
 end
 
@@ -37,7 +24,6 @@ Bidiagonal(A::AbstractMatrix, isupper::Bool)=Bidiagonal(diag(A), diag(A, isupper
 #Converting from Bidiagonal to dense Matrix
 full{T}(M::Bidiagonal{T}) = convert(Matrix{T}, M)
 convert{T}(::Type{Matrix{T}}, A::Bidiagonal{T})=diagm(A.dv) + diagm(A.ev, A.isupper?1:-1)
-promote_rule{T}(::Type{Matrix{T}}, ::Type{Bidiagonal{T}})=Matrix{T}
 promote_rule{T,S}(::Type{Matrix{T}}, ::Type{Bidiagonal{S}})=Matrix{promote_type(T,S)}
 
 #Converting from Bidiagonal to Tridiagonal
@@ -46,9 +32,19 @@ function convert{T}(::Type{Tridiagonal{T}}, A::Bidiagonal{T})
     z = zeros(T, size(A)[1]-1)
     A.isupper ? Tridiagonal(A.ev, A.dv, z) : Tridiagonal(z, A.dv, A.ev)
 end
-promote_rule{T}(::Type{Tridiagonal{T}}, ::Type{Bidiagonal{T}})=Tridiagonal{T}
 promote_rule{T,S}(::Type{Tridiagonal{T}}, ::Type{Bidiagonal{S}})=Tridiagonal{promote_type(T,S)}
 
+#################
+# BLAS routines #
+#################
+
+#Singular values
+svdvals{T<:BlasReal}(M::Bidiagonal{T})=LAPACK.bdsdc!(M.isupper?'U':'L', 'N', copy(M.dv), copy(M.ev))
+svd    {T<:BlasReal}(M::Bidiagonal{T})=LAPACK.bdsdc!(M.isupper?'U':'L', 'I', copy(M.dv), copy(M.ev))
+
+####################
+# Generic routines #
+####################
 
 function show(io::IO, M::Bidiagonal)
     println(io, summary(M), ":")
@@ -62,14 +58,30 @@ size(M::Bidiagonal) = (length(M.dv), length(M.dv))
 size(M::Bidiagonal, d::Integer) = d<1 ? error("dimension out of range") : (d<=2 ? length(M.dv) : 1)
 
 #Elementary operations
-copy(M::Bidiagonal) = Bidiagonal(copy(M.dv), copy(M.ev), copy(M.isupper))
-round(M::Bidiagonal) = Bidiagonal(round(M.dv), round(M.ev), M.isupper)
-iround(M::Bidiagonal) = Bidiagonal(iround(M.dv), iround(M.ev), M.isupper)
+for func in (conj, copy, round, iround)
+    func(M::Bidiagonal) = Bidiagonal(func(M.dv), func(M.ev), M.isupper)
+end
 
-conj(M::Bidiagonal) = Bidiagonal(conj(M.dv), conj(M.ev), M.isupper)
 transpose(M::Bidiagonal) = Bidiagonal(M.dv, M.ev, !M.isupper)
 ctranspose(M::Bidiagonal) = Bidiagonal(conj(M.dv), conj(M.ev), !M.isupper)
 
+istriu(M::Bidiagonal) = M.isupper
+istril(M::Bidiagonal) = !M.isupper
+
+function diag{T}(M::Bidiagonal{T}, n::Integer=0)
+    if n==0
+        return M.dv
+    elseif n==1
+        return M.isupper ? M.ev : zeros(T, size(M,1)-1)
+    elseif n==-1
+        return M.isupper ? zeros(T, size(M,1)-1) : M.ev
+    elseif -size(M,1)<n<size(M,1)
+        return zeros(T, size(M,1)-abs(n))
+    else
+        throw(BoundsError())
+    end
+end
+
 function +(A::Bidiagonal, B::Bidiagonal)
     if A.isupper==B.isupper
         Bidiagonal(A.dv+B.dv, A.ev+B.ev, A.isupper)
@@ -92,17 +104,101 @@ end
 /(A::Bidiagonal, B::Number) = Bidiagonal(A.dv/B, A.ev/B, A.isupper)
 ==(A::Bidiagonal, B::Bidiagonal) = (A.dv==B.dv) && (A.ev==B.ev) && (A.isupper==B.isupper)
 
-SpecialMatrix = Union(Diagonal, Bidiagonal, SymTridiagonal, Tridiagonal)
+SpecialMatrix = Union(Diagonal, Bidiagonal, SymTridiagonal, Tridiagonal, Triangular)
 *(A::SpecialMatrix, B::SpecialMatrix)=full(A)*full(B)
 
-# solver uses tridiagonal gtsv! 
-function \{T<:BlasFloat}(M::Bidiagonal{T}, rhs::StridedVecOrMat{T})
-    stride(rhs, 1)==1 || solve(M, rhs) #generic fallback
-    z = zeros(size(M, 1) - 1)
-    apply(LAPACK.gtsv!, M.isupper ? (z, copy(M.dv), copy(M.ev), copy(rhs)) : (copy(M.ev), copy(M.dv), z, copy(rhs)))
+#Generic multiplication
+for func in (:*, :Ac_mul_B, :A_mul_Bc, :/, :A_rdiv_Bc)
+    @eval begin
+        ($func){T}(A::Bidiagonal{T}, B::AbstractVector{T}) = ($func)(full(A), B)
+        #($func){T}(A::AbstractArray{T}, B::Triangular{T}) = ($func)(full(A), B)
+    end
+end
+
+
+#Linear solvers
+\{T<:Number}(A::Union(Bidiagonal{T},Triangular{T}), b::AbstractVector{T}) = naivesub!(A, b, similar(b))
+\{T<:Number}(A::Union(Bidiagonal{T},Triangular{T}), B::AbstractMatrix{T}) = hcat([naivesub!(A, B[:,i], similar(B[:,i])) for i=1:size(B,2)]...)
+A_ldiv_B!(A::Union(Bidiagonal,Triangular), b::AbstractVector)=naivesub!(A,b)
+At_ldiv_B!(A::Union(Bidiagonal,Triangular), b::AbstractVector)=naivesub!(transpose(A),b)
+Ac_ldiv_B!(A::Union(Bidiagonal,Triangular), b::AbstractVector)=naivesub!(ctranspose(A),b)
+for func in (:A_ldiv_B!, :Ac_ldiv_B!, :At_ldiv_B!) @eval begin
+    function ($func)(A::Bidiagonal, B::AbstractMatrix)
+        for i=1:size(B,2)
+            ($func)(A, B[:,i])
+        end
+        B
+    end
+end end
+for func in (:A_ldiv_Bt!, :Ac_ldiv_Bt!, :At_ldiv_Bt!) @eval begin
+    function ($func)(A::Bidiagonal, B::AbstractMatrix)
+        for i=1:size(B,1)
+            ($func)(A, B[i,:])
+        end
+        B
+    end
+end end
+
+function inv(A::Union(Bidiagonal, Triangular))
+    B = eye(eltype(A), size(A, 1))
+    for i=1:size(B,2)
+        naivesub!(A, B[:,i])
+    end
+    B
+end
+
+#Generic solver using naive substitution
+function naivesub!{T}(A::Bidiagonal{T}, b::AbstractVector, x::AbstractVector=b)
+    N = size(A, 2)
+    N==length(b)==length(x) || throw(DimensionMismatch(""))
+
+    if !A.isupper #do forward substitution
+        for j = 1:N
+            x[j] = b[j]
+            j>1 && (x[j] -= A.ev[j-1] * x[j-1])
+            x[j]/= A.dv[j]==zero(T) ? throw(SingularException(j)) : A.dv[j]
+        end
+    else #do backward substitution
+        for j = N:-1:1
+            x[j] = b[j]
+            j<N && (x[j] -= A.ev[j] * x[j+1])
+            x[j]/= A.dv[j]==zero(T) ? throw(SingularException(j)) : A.dv[j]
+        end
+    end
+    x
 end
 
-#Wrap bdsdc to compute singular values and vectors
-svdvals{T<:Real}(M::Bidiagonal{T})=LAPACK.bdsdc!(M.isupper?'U':'L', 'N', copy(M.dv), copy(M.ev))
-svd    {T<:Real}(M::Bidiagonal{T})=LAPACK.bdsdc!(M.isupper?'U':'L', 'I', copy(M.dv), copy(M.ev))
+# Eigensystems
+eigvals(M::Bidiagonal) = M.dv
+function eigvecs{T}(M::Bidiagonal{T})
+    n = length(M.dv)
+    Q=Array(T, n, n)
+    blks = [0; find(x->x==0, M.ev); n]
+    if M.isupper
+        for idx_block=1:length(blks)-1, i=blks[idx_block]+1:blks[idx_block+1] #index of eigenvector
+            v=zeros(T, n)
+            v[blks[idx_block]+1] = one(T)
+            for j=blks[idx_block]+1:i-1 #Starting from j=i, eigenvector elements will be 0
+                v[j+1] = (M.dv[i]-M.dv[j])/M.ev[j] * v[j]
+            end
+            Q[:,i] = v/norm(v)
+        end
+    else
+        for idx_block=1:length(blks)-1, i=blks[idx_block+1]:-1:blks[idx_block]+1 #index of eigenvector
+            v=zeros(T, n)
+            v[blks[idx_block+1]] = one(T)
+            for j=(blks[idx_block+1]-1):-1:max(1,(i-1)) #Starting from j=i, eigenvector elements will be 0
+                v[j] = (M.dv[i]-M.dv[j+1])/M.ev[j] * v[j+1]
+            end
+            Q[:,i] = v/norm(v)
+        end
+    end
+    Q #Actually Triangular
+end
+eigfact(M::Bidiagonal) = Eigen(eigvals(M), eigvecs(M))
 
+#Singular values
+function svdfact(M::Bidiagonal, thin::Bool=true)
+    U, S, V = svd(M)
+    SVD(U, S, V')
+end

base/linalg/triangular.jl

@@ -102,25 +102,6 @@ for func in (:*, :Ac_mul_B, :A_mul_Bc, :/, :A_rdiv_Bc)
     end
 end
 
-A_ldiv_B!(A::Triangular, b::AbstractVector)=naivesub!(A,b)
-At_ldiv_B!(A::Triangular, b::AbstractVector)=naivesub!(transpose(A),b)
-Ac_ldiv_B!(A::Triangular, b::AbstractVector)=naivesub!(ctranspose(A),b)
-for func in (:A_ldiv_B!, :Ac_ldiv_B!, :At_ldiv_B!) @eval begin
-    function ($func)(A::Triangular, B::AbstractMatrix)
-        for i=1:size(B,2)
-            ($func)(A, B[:,i])
-        end
-        B
-    end
-end end
-for func in (:A_ldiv_Bt!, :Ac_ldiv_Bt!, :At_ldiv_Bt!) @eval begin
-    function ($func)(A::Triangular, B::AbstractMatrix)
-        for i=1:size(B,1)
-            ($func)(A, B[i,:])
-        end
-        B
-    end
-end end
 #Generic solver using naive substitution
 function naivesub!(A::Triangular, b::AbstractVector, x::AbstractVector=b)
     N = size(A, 2)
@@ -152,17 +133,6 @@ function naivesub!(A::Triangular, b::AbstractVector, x::AbstractVector=b)
     x
 end
 
-\{T<:Number}(A::Triangular{T}, b::AbstractVector{T}) = naivesub!(A, b, similar(b))
-\{T<:Number}(A::Triangular{T}, B::AbstractMatrix{T}) = hcat([naivesub!(A, B[:,i], similar(B[:,i])) for i=1:size(B,2)]...)
-
-function inv(A::Triangular)
-    B = eye(eltype(A), size(A, 1))
-    for i=1:size(B,2)
-        naivesub!(A, B[:,i])
-    end
-    B
-end
-
 #Generic eigensystems
 eigvals(A::Triangular) = diag(A.UL)
 det(A::Triangular) = prod(eigvals(A))

base/mpfr.jl

@@ -75,7 +75,7 @@ BigFloat(x::Integer) = BigFloat(BigInt(x))
 BigFloat(x::Union(Bool,Int8,Int16,Int32)) = BigFloat(convert(Clong,x))
 BigFloat(x::Union(Uint8,Uint16,Uint32)) = BigFloat(convert(Culong,x))
 
-BigFloat(x::Float32) = BigFloat(float64(x))
+BigFloat(x::Union(Float16,Float32)) = BigFloat(float64(x))
 BigFloat(x::Rational) = BigFloat(num(x)) / BigFloat(den(x))
 
 convert(::Type{Rational}, x::BigFloat) = convert(Rational{BigInt}, x)