Skip to content

Commit

Permalink
Merge pull request #5277 from JuliaLang/cjh/solve-bidiag
Browse files Browse the repository at this point in the history 
Provides specialized Bidiagonal solvers
  • Loading branch information
@jiahao
jiahao committed Jan 9, 2014
2 parents 7c2301f + cc89594 commit 891224f
Show file tree
Hide file tree
Showing 5 changed files with 181 additions and 104 deletions.
5 changes: 4 additions & 1 deletion NEWS.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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])

Expand Down Expand Up @@ -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
==========================
Expand Down
170 changes: 133 additions & 37 deletions base/linalg/bidiag.jl
View file
Original file line number Diff line number Diff line change
@@ -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

Expand All @@ -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
Expand All @@ -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), ":")
Expand All @@ -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)
Expand All @@ -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
30 changes: 0 additions & 30 deletions base/linalg/triangular.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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))
Expand Down
2 changes: 1 addition & 1 deletion base/mpfr.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
Loading

0 comments on commit 891224f

Please sign in to comment.