Skip to content

Commit

Permalink
constructors for bidiagonal, diagonal, tridiagonal for consistency; a…
Browse files Browse the repository at this point in the history 
…lso accept Range inputs
  • Loading branch information
@artkuo
artkuo committed Nov 24, 2015
1 parent 6472a57 commit e0d75d7
Show file tree
Hide file tree
Showing 6 changed files with 35 additions and 5 deletions.
6 changes: 3 additions & 3 deletions base/linalg/bidiag.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,8 +12,8 @@ type Bidiagonal{T} <: AbstractMatrix{T}
new(dv, ev, isupper)
end
end
Bidiagonal{T}(dv::AbstractVector{T}, ev::AbstractVector{T}, isupper::Bool) = Bidiagonal{T}(dv, ev, isupper)
Bidiagonal{T}(dv::AbstractVector{T}, ev::AbstractVector{T}) = throw(ArgumentError("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}(collect(dv), collect(ev), isupper)

This comment has been minimized.

Sign in to view

Copy link
@fredrikekre

fredrikekre Jul 7, 2017

Member

This makes copies for Vector input, was that intentional? That is not the case for Diagonal for instance:

julia> a = rand(2); b = rand(1);

julia> D = Diagonal(a); B = Bidiagonal(a, b, true); fill!(a,0);

julia> D
B
2×2 Diagonal{Float64}:
 0.0   ⋅ 
  ⋅   0.0

julia> B
2×2 Bidiagonal{Float64}:
 0.609149  0.596931 
  ⋅        0.0392308

This comment has been minimized.

Sign in to view

Copy link
@Sacha0

Sacha0 Jul 8, 2017 •
edited
Loading

Member

Chances are the copying behavior was unintentional. I imagine the collects were introduced to handle ranges?

This comment has been minimized.

Sign in to view

Copy link
@andreasnoack

andreasnoack Jul 11, 2017

Member

Yes. I think that was the reason.

This comment has been minimized.

Sign in to view

Copy link
@artkuo

artkuo Jul 11, 2017

Author Contributor

Yes, main change was to allow AbstractVectors to be used. This necessitated collect. Which was consistent with your expectations, D or B?

For best consistency, it could be argued that everything should use a copy, but that would introduce inefficiency for fairly "low-level" code. An opposing type of consistency would be for everything to use a reference to (or view of) underlying a. But that would also be hard to enforce across all kinds of "constructed" matrices. Easy for Diagonal, but perhaps not as much for other matrices, e.g. Vandermonde, Toeplitz, Hankel, etc. At some point an actual dense matrix will be constructed, and nobody would expect a reference.

My experience as a user is that I need to be wary of modifying a and its possible effects on things built from a. If I want D to be untouched, I need to use diagonal(copy(a)). My own bugs are usually from not remembering that (in both Julia and Python). I don't usually expect the other behavior, where a can affect B after the fact. Again, this is a question of what you expect.
Bidiagonal(dv::AbstractVector, ev::AbstractVector) = throw(ArgumentError("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
Bidiagonal(dv::AbstractVector, ev::AbstractVector, uplo::Char) = begin
Expand All @@ -24,7 +24,7 @@ Bidiagonal(dv::AbstractVector, ev::AbstractVector, uplo::Char) = begin
else
throw(ArgumentError("Bidiagonal uplo argument must be upper 'U' or lower 'L', got $(repr(uplo))"))
end
Bidiagonal(copy(dv), copy(ev), isupper)
Bidiagonal(collect(dv), collect(ev), isupper)
end
function Bidiagonal{Td,Te}(dv::AbstractVector{Td}, ev::AbstractVector{Te}, isupper::Bool)
T = promote_type(Td,Te)
Expand Down
1 change: 1 addition & 0 deletions base/linalg/diagonal.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,6 +6,7 @@ immutable Diagonal{T} <: AbstractMatrix{T}
diag::Vector{T}
end
Diagonal(A::AbstractMatrix) = Diagonal(diag(A))
Diagonal(V::AbstractVector) = Diagonal(collect(V))

convert{T}(::Type{Diagonal{T}}, D::Diagonal{T}) = D
convert{T}(::Type{Diagonal{T}}, D::Diagonal) = Diagonal{T}(convert(Vector{T}, D.diag))
Expand Down
19 changes: 17 additions & 2 deletions base/linalg/tridiag.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,7 +16,7 @@ end

SymTridiagonal{T}(dv::Vector{T}, ev::Vector{T}) = SymTridiagonal{T}(dv, ev)

function SymTridiagonal{Td,Te}(dv::Vector{Td}, ev::Vector{Te})
function SymTridiagonal{Td,Te}(dv::AbstractVector{Td}, ev::AbstractVector{Te})
T = promote_type(Td,Te)
SymTridiagonal(convert(Vector{T}, dv), convert(Vector{T}, ev))
end
Expand Down Expand Up @@ -271,17 +271,32 @@ immutable Tridiagonal{T} <: AbstractMatrix{T}
du::Vector{T} # sup-diagonal
du2::Vector{T} # supsup-diagonal for pivoting
end

# Basic constructor takes in three dense vectors of same type
function Tridiagonal{T}(dl::Vector{T}, d::Vector{T}, du::Vector{T})
n = length(d)
if (length(dl) != n-1 || length(du) != n-1)
throw(ArgumentError("Cannot make Tridiagonal from incompatible lengths of subdiagonal, diagonal and superdiagonal: ($(length(dl)), $(length(d)), $(length(du))"))
end
Tridiagonal(dl, d, du, zeros(T,n-2))
end
function Tridiagonal{Tl, Td, Tu}(dl::Vector{Tl}, d::Vector{Td}, du::Vector{Tu})

# Construct from diagonals of any abstract vector, any eltype
function Tridiagonal{Tl, Td, Tu}(dl::AbstractVector{Tl}, d::AbstractVector{Td}, du::AbstractVector{Tu})
Tridiagonal(map(v->convert(Vector{promote_type(Tl,Td,Tu)}, v), (dl, d, du))...)
end

# Provide a constructor Tridiagonal(A) similar to the triangulars, diagonal, symmetric
"""
Tridiagonal(A)
returns a `Tridiagonal` array based on (abstract) matrix `A`, using its lower diagonal,
diagonal, and upper diagonal.
"""
function Tridiagonal(A::AbstractMatrix)
return Tridiagonal(diag(A,-1), diag(A), diag(A,+1))
end

size(M::Tridiagonal) = (length(M.d), length(M.d))
function size(M::Tridiagonal, d::Integer)
if d < 1
Expand Down
3 changes: 3 additions & 0 deletions test/linalg/bidiag.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -184,6 +184,9 @@ let A = Bidiagonal([1,2,3], [0,0], true)
@test isdiag(A)
end

# test construct from range
@test Bidiagonal(1:3, 1:2, true) == [1 1 0; 0 2 2; 0 0 3]

#test promote_rule
A = Bidiagonal(ones(Float32,10),ones(Float32,9),true)
B = rand(Float64,10,10)
Expand Down
3 changes: 3 additions & 0 deletions test/linalg/diagonal.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -214,3 +214,6 @@ let d = randn(n), D = Diagonal(d)
@test inv(D) ≈ inv(full(D))
end
@test_throws SingularException inv(Diagonal(zeros(n)))

# allow construct from range
@test Diagonal(linspace(1,3,3)) == Diagonal([1.,2.,3.])
8 changes: 8 additions & 0 deletions test/linalg/tridiag.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -400,3 +400,11 @@ SymTridiagonal([1, 2], [0])^3 == [1 0; 0 8]
#test convert for SymTridiagonal
@test convert(SymTridiagonal{Float64},SymTridiagonal(ones(Float32,5),ones(Float32,4))) == SymTridiagonal(ones(Float64,5),ones(Float64,4))
@test convert(AbstractMatrix{Float64},SymTridiagonal(ones(Float32,5),ones(Float32,4))) == SymTridiagonal(ones(Float64,5),ones(Float64,4))

# Test constructors from matrix
@test SymTridiagonal([1 2 3; 2 5 6; 0 6 9]) == [1 2 0; 2 5 6; 0 6 9]
@test Tridiagonal([1 2 3; 4 5 6; 7 8 9]) == [1 2 0; 4 5 6; 0 8 9]

# Test constructors with range and other abstract vectors
@test SymTridiagonal(1:3, 1:2) == [1 1 0; 1 2 2; 0 2 3]
@test Tridiagonal(4:5, 1:3, 1:2) == [1 1 0; 4 2 2; 0 5 3]

0 comments on commit e0d75d7

Please sign in to comment.