diff --git a/NEWS.md b/NEWS.md index d3a762d67ae47..5bccc53fe3a4e 100644 --- a/NEWS.md +++ b/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 ========================== diff --git a/base/linalg/bidiag.jl b/base/linalg/bidiag.jl index 5737e27fc3924..3190ca4032078 100644 --- a/base/linalg/bidiag.jl +++ b/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)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] + jx==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 diff --git a/base/linalg/triangular.jl b/base/linalg/triangular.jl index 1d900834c46d8..bacca28ffcc15 100644 --- a/base/linalg/triangular.jl +++ b/base/linalg/triangular.jl @@ -104,25 +104,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) @@ -154,17 +135,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)) diff --git a/base/mpfr.jl b/base/mpfr.jl index 629250475e130..e69f6b8f95abc 100644 --- a/base/mpfr.jl +++ b/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) diff --git a/test/linalg.jl b/test/linalg.jl index 7ac7e3342c3cc..20dab6b92d619 100644 --- a/test/linalg.jl +++ b/test/linalg.jl @@ -1,5 +1,5 @@ import Base.LinAlg -import Base.LinAlg: BlasFloat +import Base.LinAlg: BlasComplex, BlasFloat, BlasReal n = 10 srand(1234321) @@ -550,15 +550,14 @@ for elty in (Float32, Float64, Complex64, Complex128) end #Test equivalence of eigenvectors/singular vectors taking into account possible phase (sign) differences -function test_approx_eq_vecs(a, b) - n = size(a)[1] - @test n==size(b)[1] - elty = typeof(a[1]) - @test elty==typeof(b[1]) +function test_approx_eq_vecs{S<:Real,T<:Real}(a::StridedVecOrMat{S}, b::StridedVecOrMat{T}, error=nothing) + n = size(a, 1) + @test n==size(b,1) && size(a,2)==size(b,2) + if error==nothing error=n^2*(eps(S)+eps(T)) end for i=1:n ev1, ev2 = a[:,i], b[:,i] deviation = min(abs(norm(ev1-ev2)),abs(norm(ev1+ev2))) - @test_approx_eq_eps deviation 0.0 n^2*eps(abs(convert(elty, 1.0))) + @test_approx_eq_eps deviation 0.0 error end end @@ -605,7 +604,6 @@ for relty in (Float16, Float32, Float64, BigFloat), elty in (relty, Complex{relt end #SymTridiagonal (symmetric tridiagonal) matrices -n=5 Ainit = randn(n) Binit = randn(n-1) for elty in (Float32, Float64) @@ -632,46 +630,56 @@ for elty in (Float32, Float64) test_approx_eq_vecs(v, evecs) end - #Bidiagonal matrices -dv = randn(n) -ev = randn(n-1) -for elty in (Float32, Float64, Complex64, Complex128) - if (elty == Complex64) - dv += im*randn(n) - ev += im*randn(n-1) +for relty in (Float16, Float32, Float64, BigFloat), elty in (relty, Complex{relty}) + dv = convert(Vector{elty}, randn(n)) + ev = convert(Vector{elty}, randn(n-1)) + b = convert(Vector{elty}, randn(n)) + if (elty <: Complex) + dv += im*convert(Vector{elty}, randn(n)) + ev += im*convert(Vector{elty}, randn(n-1)) + b += im*convert(Vector{elty}, randn(n)) end for isupper in (true, false) #Test upper and lower bidiagonal matrices - T = Bidiagonal{elty}(dv, ev, isupper) + T = Bidiagonal(dv, ev, isupper) - @test size(T, 1) == n + @test size(T, 1) == size(T, 2) == n @test size(T) == (n, n) @test full(T) == diagm(dv) + diagm(ev, isupper?1:-1) @test Bidiagonal(full(T), isupper) == T z = zeros(elty, n) # idempotent tests - @test conj(conj(T)) == T - @test transpose(transpose(T)) == T - @test ctranspose(ctranspose(T)) == T + for func in (conj, transpose, ctranspose) + @test func(func(T)) == T + end - if (elty <: Real) - Tfull = full(T) + #Linear solver + Tfull = full(T) + condT = cond(complex128(Tfull)) + x = T \ b + tx = Tfull \ b + @test norm(x-tx,Inf) <= 4*condT*max(eps()*norm(tx,Inf), eps(relty)*norm(x,Inf)) + + #Test eigenvalues/vectors + d1, v1 = eig(T) + d2, v2 = eig((elty<:Complex?complex128:float64)(Tfull)) + @test_approx_eq isupper?d1:reverse(d1) d2 + if elty <: Real + test_approx_eq_vecs(v1, isupper?v2:v2[:,n:-1:1]) + end + + if (elty <: BlasReal) #Test singular values/vectors @test_approx_eq svdvals(Tfull) svdvals(T) u1, d1, v1 = svd(Tfull) u2, d2, v2 = svd(T) @test_approx_eq d1 d2 - test_approx_eq_vecs(u1, u2) - test_approx_eq_vecs(v1, v2) - - #Test eigenvalues/vectors - #d1, v1 = eig(Tfull) - #d2, v2 = eigvals(T), eigvecs(T) - #@test_approx_eq d1 d2 - #test_approx_eq_vecs(v1, v2) - #@test_approx_eq_eps 0 norm(v1 * diagm(d1) * inv(v1) - Tfull) eps(elty)*n*(n+1) - #@test_approx_eq_eps 0 norm(v2 * diagm(d2) * inv(v2) - Tfull) eps(elty)*n*(n+1) + if elty <: Real + test_approx_eq_vecs(u1, u2) + test_approx_eq_vecs(v1, v2) + end + @test_approx_eq_eps 0 normfro(u2*diagm(d2)*v2'-Tfull) n*max(n^2*eps(relty), normfro(u1*diagm(d1)*v1'-Tfull)) end end end @@ -692,8 +700,8 @@ for relty in (Float16, Float32, Float64, BigFloat), elty in (relty, Complex{relt @test_approx_eq_eps D*v DM*v n*eps(relty)*(elty<:Complex ? 2:1) @test_approx_eq_eps D*U DM*U n^2*eps(relty)*(elty<:Complex ? 2:1) if relty != BigFloat - @test_approx_eq_eps D\v DM\v n*eps(relty)*(elty<:Complex ? 2:1) - @test_approx_eq_eps D\U DM\U n^2*eps(relty)*(elty<:Complex ? 2:1) + @test_approx_eq_eps D\v DM\v 2n^2*eps(relty)*(elty<:Complex ? 2:1) + @test_approx_eq_eps D\U DM\U 2n^3*eps(relty)*(elty<:Complex ? 2:1) end for func in (det, trace) @test_approx_eq_eps func(D) func(DM) n^2*eps(relty) @@ -703,7 +711,7 @@ for relty in (Float16, Float32, Float64, BigFloat), elty in (relty, Complex{relt @test_approx_eq_eps func(D) func(DM) n^2*eps(relty) end end - if elty <: Complex && relty <:BlasFloat + if elty <: BlasComplex for func in (logdet, sqrtm) @test_approx_eq_eps func(D) func(DM) n^2*eps(relty) end