Skip to content

Commit

Permalink
parametrize SymTridiagonal on the wrapped vector type
Browse files Browse the repository at this point in the history
  • Loading branch information
fredrikekre committed Aug 3, 2017
1 parent 42b521f commit a3d50cf
Show file tree
Hide file tree
Showing 5 changed files with 74 additions and 35 deletions.
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
`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}
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 @@ -28,6 +28,19 @@ for elty in (Float32, Float64, Complex64, Complex128, Int)
v = convert(Vector{elty}, v)
B = convert(Matrix{elty}, B)
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

ε = eps(abs2(float(one(elty))))
T = Tridiagonal(dl, d, du)
Ts = SymTridiagonal(d, dl)
Expand Down Expand Up @@ -251,7 +264,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

0 comments on commit a3d50cf

Please sign in to comment.