Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Provides Bidiagonal matrix type and specialized SVD #2713

Merged
merged 1 commit into from
Mar 31, 2013
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions base/exports.jl
Original file line number Diff line number Diff line change
Expand Up @@ -22,6 +22,7 @@ export
Array,
Associative,
AsyncStream,
Bidiagonal,
BitArray,
BigFloat,
BigInt,
Expand Down
2 changes: 2 additions & 0 deletions base/linalg.jl
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,7 @@ export
BunchKaufman,
SymTridiagonal,
Tridiagonal,
Bidiagonal,
Woodbury,
Factorization,
BunchKaufman,
Expand Down Expand Up @@ -166,6 +167,7 @@ include("linalg/triangular.jl")
include("linalg/hermitian.jl")
include("linalg/woodbury.jl")
include("linalg/tridiag.jl")
include("linalg/bidiag.jl")
include("linalg/diagonal.jl")
include("linalg/rectfullpacked.jl")

Expand Down
120 changes: 120 additions & 0 deletions base/linalg/bidiag.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,120 @@
#### Specialized matrix types ####

## Bidiagonal matrices


type Bidiagonal{T} <: AbstractMatrix{T}
dv::Vector{T} # diagonal
ev::Vector{T} # sub/super diagonal
isupper::Bool # is upper bidiagonal (true) or lower (false)
function Bidiagonal{T}(dv::Vector{T}, ev::Vector{T}, isupper::Bool)
if length(ev)!=length(dv)-1 error("dimension mismatch") end
new(dv, ev, isupper)
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.")


#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(string("Bidiagonal can only be upper 'U' or lower 'L' but you said '", uplo, "'"))
end
Bidiagonal{T}(copy(dv), copy(ev), isupper)
end

function Bidiagonal{Td<:Number,Te<:Number}(dv::Vector{Td}, ev::Vector{Te}, isupper::Bool)
T = promote(Td,Te)
Bidiagonal(convert(Vector{T}, dv), convert(Vector{T}, ev), isupper)
end

Bidiagonal(A::AbstractMatrix, isupper::Bool)=Bidiagonal(diag(A), diag(A, isupper?1:-1), 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
Tridiagonal{T}(M::Bidiagonal{T}) = convert(Tridiagonal{T}, M)
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)}


function show(io::IO, M::Bidiagonal)
println(io, summary(M), ":")
print(io, "diag: ")
print_matrix(io, (M.dv)')
print(io, M.isupper?"\n sup: ":"\n sub: ")
print_matrix(io, (M.ev)')
end

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)

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)

function +(A::Bidiagonal, B::Bidiagonal)
if A.isupper==B.isupper
Bidiagonal(A.dv+B.dv, A.ev+B.ev, A.isupper)
else #return tridiagonal
if A.isupper #&& !B.isupper
Tridiagonal(B.ev,A.dv+B.dv,A.ev)
else
Tridiagonal(A.ev,A.dv+B.dv,B.ev)
end
end
end

function -(A::Bidiagonal, B::Bidiagonal)
if A.isupper==B.isupper
Bidiagonal(A.dv-B.dv, A.ev-B.ev, A.isupper)
else #return tridiagonal
if A.isupper #&& !B.isupper
Tridiagonal(-B.ev,A.dv-B.dv,A.ev)
else
Tridiagonal(A.ev,A.dv-B.dv,-B.ev)
end
end
end

-(A::Bidiagonal)=Bidiagonal(-A.dv,-A.ev)
#XXX Returns dense matrix but really should be banded
*(A::Bidiagonal, B::Bidiagonal) = full(A)*full(B)
==(A::Bidiagonal, B::Bidiagonal) = (A.dv==B.dv) && (A.ev==B.ev) && (A.isupper==B.isupper)

# solver uses tridiagonal gtsv!
function \{T<:BlasFloat}(M::Bidiagonal{T}, rhs::StridedVecOrMat{T})
if stride(rhs, 1) == 1
z = zeros(size(M)[1])
if M.isupper
return LAPACK.gtsv!(z, copy(M.dv), copy(M.ev), copy(rhs))
else
return LAPACK.gtsv!(copy(M.ev), copy(M.dv), z, copy(rhs))
end
end
solve(M, rhs) # use the Julia "fallback"
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))

109 changes: 109 additions & 0 deletions base/linalg/lapack.jl
Original file line number Diff line number Diff line change
Expand Up @@ -1659,6 +1659,115 @@ for (syconv, syev, sysv, sytrf, sytri, sytrs, elty, relty) in
end
end



#Find the leading dimension
ld = x->max(1,stride(x,2))
function validate(uplo)
if !(uplo=='U' || uplo=='L')
error(string("Invalid UPLO: must be 'U' or 'L' but you said", uplo))
end
end
## (BD) Bidiagonal matrices - singular value decomposition
for (bdsqr, relty, elty) in
((:dbdsqr_,:Float64,:Float64),
(:sbdsqr_,:Float32,:Float32),
(:zbdsqr_,:Float64,:Complex128),
(:cbdsqr_,:Float32,:Complex64))
@eval begin
#*> DBDSQR computes the singular values and, optionally, the right and/or
#*> left singular vectors from the singular value decomposition (SVD) of
#*> a real N-by-N (upper or lower) bidiagonal matrix B using the implicit
#*> zero-shift QR algorithm.
function bdsqr!(uplo::BlasChar, d::Vector{$relty}, e_::Vector{$relty},
vt::StridedMatrix{$elty}, u::StridedMatrix{$elty}, c::StridedMatrix{$elty})

validate(uplo)
n = length(d)
if length(e_) != n-1 throw(DimensionMismatch("bdsqr!")) end
ncvt, nru, ncc = size(vt)[2], size(u)[1], size(c)[2]
ldvt, ldu, ldc = ld(vt), ld(u), ld(c)
work=Array($elty, 4n)
info=Array(BlasInt,1)

ccall(($(string(bdsqr)),liblapack), Void,
(Ptr{BlasChar}, Ptr{BlasInt}, Ptr{BlasInt}, Ptr{BlasInt}, Ptr{BlasInt},
Ptr{$elty}, Ptr{$elty}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty},
Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}),
&uplo, &n, ncvt, &nru, &ncc,
d, e_, vt, &ldvt, u,
&ldu, c, &ldc, work, info)

if info[1] != 0 throw(LAPACKException(info[1])) end
d, vt, u, c #singular values in descending order, P**T * VT, U * Q, Q**T * C
end
end
end

#Defined only for real types
for (bdsdc, elty) in
((:dbdsdc_,:Float64),
(:sbdsdc_,:Float32))
@eval begin
#* DBDSDC computes the singular value decomposition (SVD) of a real
#* N-by-N (upper or lower) bidiagonal matrix B: B = U * S * VT,
#* using a divide and conquer method
#* .. Scalar Arguments ..
# CHARACTER COMPQ, UPLO
# INTEGER INFO, LDU, LDVT, N
#* ..
#* .. Array Arguments ..
# INTEGER IQ( * ), IWORK( * )
# DOUBLE PRECISION D( * ), E( * ), Q( * ), U( LDU, * ),
# $ VT( LDVT, * ), WORK( * )
function bdsdc!(uplo::BlasChar, compq::BlasChar, d::Vector{$elty}, e_::Vector{$elty})
validate(uplo)
n, ldiq, ldq, ldu, ldvt = length(d), 1, 1, 1, 1
if compq == 'N'
lwork = 4n
elseif compq == 'P'
warn("COMPQ='P' is not tested")
#TODO turn this into an actual LAPACK call
#smlsiz=ilaenv(9, $elty==:Float64 ? 'dbdsqr' : 'sbdsqr', string(uplo, compq), n,n,n,n)
smlsiz=100 #For now, completely overkill
ldq = n*(11+2*smlsiz+8*int(log((n/(smlsiz+1)))/log(2)))
ldiq = n*(3+3*int(log(n/(smlsiz+1))/log(2)))
lwork = 6n
elseif compq == 'I'
ldvt=ldu=max(1, n)
lwork=3*n^2 + 4n
else
error(string("Invalid COMPQ. Valid choices are 'N', 'P' or 'I' but you said '",compq,"'"))
end
u = Array($elty, (ldu, n))
vt= Array($elty, (ldvt, n))
q = Array($elty, ldq)
iq= Array(BlasInt, ldiq)
work =Array($elty, lwork)
iwork=Array(BlasInt, 7n)
info =Array(BlasInt, 1)
ccall(($(string(bdsdc)),liblapack), Void,
(Ptr{BlasChar}, Ptr{BlasChar}, Ptr{BlasInt}, Ptr{$elty}, Ptr{$elty},
Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt},
Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{BlasInt}),
&uplo, &compq, &n, d, e_,
u, &ldu, vt, &ldvt,
q, iq, work, iwork, info)

if info[1] != 0 throw(LAPACKException(info[1])) end
if compq == 'N'
d
elseif compq == 'P'
d, q, iq
else #compq == 'I'
u, d, vt'
end
end
end
end



# New symmetric eigen solver
for (syevr, elty) in
((:dsyevr_,:Float64),
Expand Down
18 changes: 5 additions & 13 deletions base/linalg/tridiag.jl
Original file line number Diff line number Diff line change
Expand Up @@ -23,17 +23,9 @@ function SymTridiagonal{Td<:Number,Te<:Number}(dv::Vector{Td}, ev::Vector{Te})
SymTridiagonal(convert(Vector{T}, dv), convert(Vector{T}, ev))
end

SymTridiagonal(A::AbstractMatrix) = SymTridiagonal(diag(A), diag(A,1))

SymTridiagonal(M::AbstractMatrix) = diag(A,1)==diag(A,-1)?SymTridiagonal(diag(A), diag(A,1)):error("Matrix is not symmetric, cannot convert to SymTridiagonal")
full{T}(M::SymTridiagonal{T}) = convert(Matrix{T}, M)
function convert{T}(::Type{Matrix{T}}, S::SymTridiagonal{T})
M = diagm(S.dv)
for i in 1:length(S.ev)
j = i + 1
M[i,j] = M[j,i] = S.ev[i]
end
M
end
convert{T}(::Type{Matrix{T}}, M::SymTridiagonal{T})=diagm(M.dv)+diagm(M.ev,-1)+diagm(M.ev,1)

function show(io::IO, S::SymTridiagonal)
println(io, summary(S), ":")
Expand Down Expand Up @@ -96,8 +88,8 @@ end

function Tridiagonal{T<:Number}(dl::Vector{T}, d::Vector{T}, du::Vector{T})
N = length(d)
if length(dl) != N-1 || length(du) != N-1
error("The sub- and super-diagonals must have length N-1")
if (length(dl) != N-1 || length(du) != N-1)
error(string("Cannot make Tridiagonal from incompatible lengths of subdiagonal, diagonal and superdiagonal: (", length(dl), ", ", length(d), ", ", length(du),")"))
end
M = Tridiagonal{T}(N)
M.dl = copy(dl)
Expand Down Expand Up @@ -160,7 +152,7 @@ ctranspose(M::Tridiagonal) = conj(transpose(M))
==(A::SymTridiagonal, B::SymTridiagonal) = B==A

# Elementary operations that mix Tridiagonal and SymTridiagonal matrices
Tridiagonal(A::SymTridiagonal) = Tridiagonal(A.dv, A.ev, A.dv)
convert(::Type{Tridiagonal}, A::SymTridiagonal) = Tridiagonal(A.ev, A.dv, A.ev)
+(A::Tridiagonal, B::SymTridiagonal) = Tridiagonal(A.dl+B.ev, A.d+B.dv, A.du+B.ev)
+(A::SymTridiagonal, B::Tridiagonal) = Tridiagonal(A.ev+B.dl, A.dv+B.d, A.ev+B.du)
-(A::Tridiagonal, B::SymTridiagonal) = Tridiagonal(A.dl-B.ev, A.d-B.dv, A.du-B.ev)
Expand Down
7 changes: 7 additions & 0 deletions doc/helpdb.jl
Original file line number Diff line number Diff line change
Expand Up @@ -5094,6 +5094,13 @@
"),

("Linear Algebra","","Bidiagonal","Bidiagonal(dv, ev, isupper)
Construct an upper (isupper=true) or lower (isupper=false) bidiagonal matrix
from the given diagonal (dv) and off-diagonal (ev) vectors.
"),

("Linear Algebra","","Woodbury","Woodbury(A, U, C, V)
Construct a matrix in a form suitable for applying the Woodbury
Expand Down
5 changes: 5 additions & 0 deletions doc/stdlib/linalg.rst
Original file line number Diff line number Diff line change
Expand Up @@ -175,6 +175,11 @@ Linear algebra functions in Julia are largely implemented by calling functions f

Construct a tridiagonal matrix from the lower diagonal, diagonal, and upper diagonal

.. function:: Bidiagonal(dv, ev, isupper)

Constructs an upper (isupper=true) or lower (isupper=false) bidiagonal matrix
using the given diagonal (dv) and off-diagonal (ev) vectors

.. function:: Woodbury(A, U, C, V)

Construct a matrix in a form suitable for applying the Woodbury matrix identity
Expand Down
Loading