Merge pull request #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