diff --git a/base/exports.jl b/base/exports.jl index 5a646bdf255c9..09ea9466190cf 100644 --- a/base/exports.jl +++ b/base/exports.jl @@ -22,6 +22,7 @@ export Array, Associative, AsyncStream, + Bidiagonal, BitArray, BigFloat, BigInt, diff --git a/base/linalg.jl b/base/linalg.jl index ede3230a5f00c..974e4c1b33b6b 100644 --- a/base/linalg.jl +++ b/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") diff --git a/base/linalg/bidiag.jl b/base/linalg/bidiag.jl new file mode 100644 index 0000000000000..6117b114e58ec --- /dev/null +++ b/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)) + diff --git a/base/linalg/lapack.jl b/base/linalg/lapack.jl index 491f911d1344a..4efb853757323 100644 --- a/base/linalg/lapack.jl +++ b/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), diff --git a/base/linalg/tridiag.jl b/base/linalg/tridiag.jl index 4df1af7b18834..c79ee7e588e0a 100644 --- a/base/linalg/tridiag.jl +++ b/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) diff --git a/doc/helpdb.jl b/doc/helpdb.jl index a6b5974b1d0ad..09a6cf8c27d6f 100644 --- a/doc/helpdb.jl +++ b/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 diff --git a/doc/stdlib/linalg.rst b/doc/stdlib/linalg.rst index 65e225ebd9eed..bc9799ddb7285 100644 --- a/doc/stdlib/linalg.rst +++ b/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 diff --git a/test/linalg.jl b/test/linalg.jl index 48d6fc34ce363..378d72cf0bcf3 100644 --- a/test/linalg.jl +++ b/test/linalg.jl @@ -383,6 +383,19 @@ for elty in (Float32, Float64, Complex64, Complex128) @test_approx_eq vals LinAlg.LAPACK.syev!('N','U',copy(Asym)) 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]) + @assert elty==typeof(b[1]) + 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))) + end +end + #LAPACK tests for symmetric tridiagonal matrices n=5 Ainit = randn(n) @@ -414,16 +427,47 @@ for elty in (Float32, Float64) (evecs, ifail, info)=LinAlg.LAPACK.stein!(A, B, w, iblock, isplit) @test info==0 @test all(ifail .== 0) - for i=1:n - ev1 = v[:,i] - ev2 = evecs[:,i] - deviation = min(abs(norm(ev1-ev2)),abs(norm(ev1+ev2))) - @test_approx_eq_eps deviation 0.0 n*eps(abs(convert(elty, 1.0))) - end + test_approx_eq_vecs(v, evecs) +end + +#Test bidiagonal matrices and their SVDs +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) + end + for isupper in (true, false) #Test upper and lower bidiagonal matrices + T = Bidiagonal{elty}(dv, ev, isupper) + + @test size(T, 1) == 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 + + if (elty <: Real) + #XXX If I run either of these tests separately, by themselves, things are OK. + # Enabling BOTH tests results in segfault. + # Where is the memory corruption??? + #@test_approx_eq svdvals(full(T)) svdvals(T) + u1, d1, v1 = svd(full(T)) + u2, d2, v2 = svd(T) + @test_approx_eq d1 d2 + test_approx_eq_vecs(u1, u2) + test_approx_eq_vecs(v1, v2) + end + end end + ## Issue related tests # issue 1447 let