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

parametrize SymTridiagonal on the wrapped vector type #23035

Merged
merged 3 commits into from
Aug 5, 2017
Merged
Show file tree
Hide file tree
Changes from 1 commit
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
Next Next commit
parametrize SymTridiagonal on the wrapped vector type
  • Loading branch information
fredrikekre committed Aug 4, 2017
commit 538904c79827e8649d5a59076032342b8c91a9ca
17 changes: 9 additions & 8 deletions NEWS.md
Original file line number Diff line number Diff line change
Expand Up @@ -83,9 +83,10 @@ This section lists changes that do not have deprecation warnings.
longer present. Use `first(R)` and `last(R)` to obtain
start/stop. ([#20974])

* The `Diagonal` and `Bidiagonal` type definitions have changed from `Diagonal{T}` and
`Bidiagonal{T}` to `Diagonal{T,V<:AbstractVector{T}}` and
`Bidiagonal{T,V<:AbstractVector{T}}` respectively ([#22718], [#22925]).
* The `Diagonal`, `Bidiagonal` and `SymTridiagonal` type definitions have changed from
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

assuming that merges first? I'm still not so sure about that densifying constructor over there

`Diagonal{T}`, `Bidiagonal{T}` and `SymTridiagonal{T}` to `Diagonal{T,V<:AbstractVector{T}}`,
`Bidiagonal{T,V<:AbstractVector{T}}` and `SymTridiagonal{T,V<:AbstractVector{T}}`
respectively ([#22718], [#22925], [#23035]).

* Spaces are no longer allowed between `@` and the name of a macro in a macro call ([#22868]).

Expand Down Expand Up @@ -153,9 +154,9 @@ Library improvements

* `Char`s can now be concatenated with `String`s and/or other `Char`s using `*` ([#22532]).

* `Diagonal` and `Bidiagonal` are now parameterized on the type of the wrapped vectors,
allowing `Diagonal` and `Bidiagonal` matrices with arbitrary
`AbstractVector`s ([#22718], [#22925]).
* `Diagonal`, `Bidiagonal` and `SymTridiagonal` are now parameterized on the type of the
wrapped vectors, allowing `Diagonal`, `Bidiagonal` and `SymTridiagonal` matrices with
arbitrary `AbstractVector`s ([#22718], [#22925], [#23035]).

* Mutating versions of `randperm` and `randcycle` have been added:
`randperm!` and `randcycle!` ([#22723]).
Expand Down Expand Up @@ -218,8 +219,8 @@ Deprecated or removed
* `Bidiagonal` constructors now use a `Symbol` (`:U` or `:L`) for the upper/lower
argument, instead of a `Bool` or a `Char` ([#22703]).

* `Bidiagonal` constructors that automatically converted the input vectors to
the same type are deprecated in favor of explicit conversion ([#22925]).
* `Bidiagonal` and `SymTridiagonal` constructors that automatically converted the input
vectors to the same type are deprecated in favor of explicit conversion ([#22925], [#23035]).

* Calling `nfields` on a type to find out how many fields its instances have is deprecated.
Use `fieldcount` instead. Use `nfields` only to get the number of fields in a specific object ([#22350]).
Expand Down
9 changes: 9 additions & 0 deletions base/deprecated.jl
Original file line number Diff line number Diff line change
Expand Up @@ -1604,6 +1604,15 @@ function Bidiagonal(dv::AbstractVector{T}, ev::AbstractVector{S}, uplo::Symbol)
Bidiagonal(convert(Vector{R}, dv), convert(Vector{R}, ev), uplo)
end

# PR #23035
# also uncomment constructor tests in test/linalg/tridiag.jl
function SymTridiagonal(dv::AbstractVector{T}, ev::AbstractVector{S}) where {T,S}
Copy link
Contributor

@tkelman tkelman Aug 1, 2017

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

could also use a deprecation bullet in news (or could add it for both in whichever merges second)

depwarn(string("SymTridiagonal(dv::AbstractVector{T}, ev::AbstractVector{S}) ",
"where {T,S} is deprecated; convert both vectors to the same type instead."), :SymTridiagonal)
R = promote_type(T, S)
SymTridiagonal(convert(Vector{R}, dv), convert(Vector{R}, ev))
end

# END 0.7 deprecations

# BEGIN 1.0 deprecations
Expand Down
66 changes: 41 additions & 25 deletions base/linalg/tridiag.jl
Original file line number Diff line number Diff line change
Expand Up @@ -3,68 +3,84 @@
#### Specialized matrix types ####

## (complex) symmetric tridiagonal matrices
struct SymTridiagonal{T} <: AbstractMatrix{T}
dv::Vector{T} # diagonal
ev::Vector{T} # subdiagonal
function SymTridiagonal{T}(dv::Vector{T}, ev::Vector{T}) where T
struct SymTridiagonal{T,V<:AbstractVector{T}} <: AbstractMatrix{T}
dv::V # diagonal
ev::V # subdiagonal
function SymTridiagonal{T}(dv::V, ev::V) where {T,V<:AbstractVector{T}}
if !(length(dv) - 1 <= length(ev) <= length(dv))
throw(DimensionMismatch("subdiagonal has wrong length. Has length $(length(ev)), but should be either $(length(dv) - 1) or $(length(dv))."))
end
new(dv,ev)
new{T,V}(dv,ev)
end
end

"""
SymTridiagonal(dv, ev)
SymTridiagonal(dv::V, ev::V) where V <: AbstractVector

Construct a symmetric tridiagonal matrix from the diagonal and first sub/super-diagonal,
respectively. The result is of type `SymTridiagonal` and provides efficient specialized
eigensolvers, but may be converted into a regular matrix with
[`convert(Array, _)`](@ref) (or `Array(_)` for short).
Construct a symmetric tridiagonal matrix from the diagonal (`dv`) and first
sub/super-diagonal (`ev`), respectively. The result is of type `SymTridiagonal`
and provides efficient specialized eigensolvers, but may be converted into a
regular matrix with [`convert(Array, _)`](@ref) (or `Array(_)` for short).

# Examples
```jldoctest
julia> dv = [1; 2; 3; 4]
julia> dv = [1, 2, 3, 4]
4-element Array{Int64,1}:
1
2
3
4

julia> ev = [7; 8; 9]
julia> ev = [7, 8, 9]
3-element Array{Int64,1}:
7
8
9

julia> SymTridiagonal(dv, ev)
4×4 SymTridiagonal{Int64}:
4×4 SymTridiagonal{Int64,Array{Int64,1}}:
1 7 ⋅ ⋅
7 2 8 ⋅
⋅ 8 3 9
⋅ ⋅ 9 4
```
"""
SymTridiagonal(dv::Vector{T}, ev::Vector{T}) where {T} = SymTridiagonal{T}(dv, ev)
SymTridiagonal(dv::V, ev::V) where {T,V<:AbstractVector{T}} = SymTridiagonal{T}(dv, ev)

function SymTridiagonal(dv::AbstractVector{Td}, ev::AbstractVector{Te}) where {Td,Te}
T = promote_type(Td,Te)
SymTridiagonal(convert(Vector{T}, dv), convert(Vector{T}, ev))
end
"""
SymTridiagonal(A::AbstractMatrix)

Construct a symmetric tridiagonal matrix from the diagonal and
first sub/super-diagonal, of the symmetric matrix `A`.

# Examples
```jldoctest
julia> A = [1 2 3; 2 4 5; 3 5 6]
3×3 Array{Int64,2}:
1 2 3
2 4 5
3 5 6

julia> SymTridiagonal(A)
3×3 SymTridiagonal{Int64,Array{Int64,1}}:
1 2 ⋅
2 4 5
⋅ 5 6
```
"""
function SymTridiagonal(A::AbstractMatrix)
if diag(A,1) == diag(A,-1)
SymTridiagonal(diag(A), diag(A,1))
SymTridiagonal(diag(A,0), diag(A,1))
else
throw(ArgumentError("matrix is not symmetric; cannot convert to SymTridiagonal"))
end
end

convert(::Type{SymTridiagonal{T}}, S::SymTridiagonal) where {T} =
SymTridiagonal(convert(Vector{T}, S.dv), convert(Vector{T}, S.ev))
SymTridiagonal(convert(AbstractVector{T}, S.dv), convert(AbstractVector{T}, S.ev))
convert(::Type{AbstractMatrix{T}}, S::SymTridiagonal) where {T} =
SymTridiagonal(convert(Vector{T}, S.dv), convert(Vector{T}, S.ev))
function convert(::Type{Matrix{T}}, M::SymTridiagonal{T}) where T
SymTridiagonal(convert(AbstractVector{T}, S.dv), convert(AbstractVector{T}, S.ev))
function convert(::Type{Matrix{T}}, M::SymTridiagonal) where T
n = size(M, 1)
Mf = zeros(T, n, n)
@inbounds begin
Expand Down Expand Up @@ -311,7 +327,7 @@ end
# R. Usmani, "Inversion of a tridiagonal Jacobi matrix",
# Linear Algebra and its Applications 212-213 (1994), pp.413-414
# doi:10.1016/0024-3795(94)90414-6
function inv_usmani(a::Vector{T}, b::Vector{T}, c::Vector{T}) where T
function inv_usmani(a::V, b::V, c::V) where {T,V<:AbstractVector{T}}
n = length(b)
θ = ZeroOffsetVector(zeros(T, n+1)) #principal minors of A
θ[0] = 1
Expand Down Expand Up @@ -341,7 +357,7 @@ end

#Implements the determinant using principal minors
#Inputs and reference are as above for inv_usmani()
function det_usmani(a::Vector{T}, b::Vector{T}, c::Vector{T}) where T
function det_usmani(a::V, b::V, c::V) where {T,V<:AbstractVector{T}}
n = length(b)
θa = one(T)
if n == 0
Expand Down Expand Up @@ -635,7 +651,7 @@ convert(::Type{AbstractMatrix{T}},M::Tridiagonal) where {T} = convert(Tridiagona
convert(::Type{Tridiagonal{T}}, M::SymTridiagonal{T}) where {T} = Tridiagonal(M)
function convert(::Type{SymTridiagonal{T}}, M::Tridiagonal) where T
if M.dl == M.du
return SymTridiagonal(convert(Vector{T},M.d), convert(Vector{T},M.dl))
return SymTridiagonal{T}(convert(AbstractVector{T},M.d), convert(AbstractVector{T},M.dl))
else
throw(ArgumentError("Tridiagonal is not symmetric, cannot convert to SymTridiagonal"))
end
Expand Down
15 changes: 14 additions & 1 deletion test/linalg/tridiag.jl
Original file line number Diff line number Diff line change
Expand Up @@ -34,6 +34,19 @@ B = randn(n,2)
F[i,i+1] = du[i]
F[i+1,i] = dl[i]
end

@testset "constructor" begin
for (x, y) in ((d, dl), (GenericArray(d), GenericArray(dl)))
ST = (SymTridiagonal(x, y))::SymTridiagonal{elty, typeof(x)}
@test ST == Matrix(ST)
@test ST.dv === x
@test ST.ev === y
end
# enable when deprecations for 0.7 are dropped
# @test_throws MethodError SymTridiagonal(dv, GenericArray(ev))
# @test_throws MethodError SymTridiagonal(GenericArray(dv), ev)
end

@testset "size and Array" begin
@test_throws ArgumentError size(Ts,0)
@test size(Ts,3) == 1
Expand Down Expand Up @@ -216,7 +229,7 @@ let n = 12 #Size of matrix problem to test
b += im*convert(Vector{elty}, randn(n-1))
end

@test_throws DimensionMismatch SymTridiagonal(a, ones(n+1))
@test_throws DimensionMismatch SymTridiagonal(a, ones(elty, n+1))
@test_throws ArgumentError SymTridiagonal(rand(n,n))

A = SymTridiagonal(a, b)
Expand Down
2 changes: 1 addition & 1 deletion test/show.jl
Original file line number Diff line number Diff line change
Expand Up @@ -558,7 +558,7 @@ A = reshape(1:16,4,4)
@test replstr(Diagonal(A)) == "4×4 Diagonal{$(Int),Array{$(Int),1}}:\n 1 ⋅ ⋅ ⋅\n ⋅ 6 ⋅ ⋅\n ⋅ ⋅ 11 ⋅\n ⋅ ⋅ ⋅ 16"
@test replstr(Bidiagonal(A,:U)) == "4×4 Bidiagonal{$(Int),Array{$(Int),1}}:\n 1 5 ⋅ ⋅\n ⋅ 6 10 ⋅\n ⋅ ⋅ 11 15\n ⋅ ⋅ ⋅ 16"
@test replstr(Bidiagonal(A,:L)) == "4×4 Bidiagonal{$(Int),Array{$(Int),1}}:\n 1 ⋅ ⋅ ⋅\n 2 6 ⋅ ⋅\n ⋅ 7 11 ⋅\n ⋅ ⋅ 12 16"
@test replstr(SymTridiagonal(A+A')) == "4×4 SymTridiagonal{$Int}:\n 2 7 ⋅ ⋅\n 7 12 17 ⋅\n ⋅ 17 22 27\n ⋅ ⋅ 27 32"
@test replstr(SymTridiagonal(A+A')) == "4×4 SymTridiagonal{$(Int),Array{$(Int),1}}:\n 2 7 ⋅ ⋅\n 7 12 17 ⋅\n ⋅ 17 22 27\n ⋅ ⋅ 27 32"
@test replstr(Tridiagonal(diag(A,-1),diag(A),diag(A,+1))) == "4×4 Tridiagonal{$Int}:\n 1 5 ⋅ ⋅\n 2 6 10 ⋅\n ⋅ 7 11 15\n ⋅ ⋅ 12 16"
@test replstr(UpperTriangular(copy(A))) == "4×4 UpperTriangular{$Int,Array{$Int,2}}:\n 1 5 9 13\n ⋅ 6 10 14\n ⋅ ⋅ 11 15\n ⋅ ⋅ ⋅ 16"
@test replstr(LowerTriangular(copy(A))) == "4×4 LowerTriangular{$Int,Array{$Int,2}}:\n 1 ⋅ ⋅ ⋅\n 2 6 ⋅ ⋅\n 3 7 11 ⋅\n 4 8 12 16"
Expand Down