Merge pull request JuliaLang#2713 from jiahao/bdsqr

Provides Bidiagonal matrix type and specialized SVD

base/exports.jl

@@ -22,6 +22,7 @@ export
  Array,
  Associative,
  AsyncStream,
+ Bidiagonal,
  BitArray,
  BigFloat,
  BigInt,

base/linalg.jl

@@ -8,6 +8,7 @@ export
  BunchKaufman,
  SymTridiagonal,
  Tridiagonal,
+ Bidiagonal,
  Woodbury,
  Factorization,
  BunchKaufman,
@@ -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")
 

base/linalg/bidiag.jl

@@ -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))
+

base/linalg/lapack.jl

@@ -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),

base/linalg/tridiag.jl

@@ -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), ":")
@@ -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)
@@ -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)

doc/helpdb.jl

@@ -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

doc/stdlib/linalg.rst

@@ -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