parametrize SymTridiagonal on the wrapped vector type (JuliaLang#23035)

* parametrize SymTridiagonal on the wrapped vector type

NEWS.md

@@ -86,9 +86,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]).
 
@@ -156,9 +157,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]).
@@ -221,8 +222,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]).

base/deprecated.jl

@@ -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

base/linalg/ldlt.jl

@@ -21,15 +21,15 @@ convert(::Type{Factorization{T}}, F::LDLt{S,U}) where {T,S,U} = convert(LDLt{T,U
 
 Same as [`ldltfact`](@ref), but saves space by overwriting the input `A`, instead of creating a copy.
 """
-function ldltfact!(S::SymTridiagonal{T}) where T<:Real
+function ldltfact!(S::SymTridiagonal{T,V}) where {T<:Real,V}
  n = size(S,1)
  d = S.dv
  e = S.ev
  @inbounds @simd for i = 1:n-1
  e[i] /= d[i]
  d[i+1] -= abs2(e[i])*d[i]
  end
- return LDLt{T,SymTridiagonal{T}}(S)
+ return LDLt{T,SymTridiagonal{T,V}}(S)
 end
 
 """
@@ -46,7 +46,7 @@ end
 
 factorize(S::SymTridiagonal) = ldltfact(S)
 
-function A_ldiv_B!(S::LDLt{T,SymTridiagonal{T}}, B::AbstractVecOrMat{T}) where T
+function A_ldiv_B!(S::LDLt{T,M}, B::AbstractVecOrMat{T}) where {T,M<:SymTridiagonal{T}}
  n, nrhs = size(B, 1), size(B, 2)
  if size(S,1) != n
  throw(DimensionMismatch("Matrix has dimensions $(size(S)) but right hand side has first dimension $n"))

base/linalg/tridiag.jl

@@ -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
@@ -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
@@ -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
@@ -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

test/linalg/lapack.jl

@@ -448,40 +448,42 @@ end
 
 @testset "ptsv" begin
  @testset for elty in (Float32, Float64, Complex64, Complex128)
- dv = real(ones(elty,10))
+ dv = ones(elty,10)
  ev = zeros(elty,9)
+ rdv = real(dv)
  A = SymTridiagonal(dv,ev)
  if elty <: Complex
  A = Tridiagonal(conj(ev),dv,ev)
  end
  B = rand(elty,10,10)
  C = copy(B)
- @test A\B ≈ LAPACK.ptsv!(dv,ev,C)
- @test_throws DimensionMismatch LAPACK.ptsv!(dv,ones(elty,10),C)
- @test_throws DimensionMismatch LAPACK.ptsv!(dv,ev,ones(elty,11,11))
+ @test A\B ≈ LAPACK.ptsv!(rdv,ev,C)
+ @test_throws DimensionMismatch LAPACK.ptsv!(rdv,ones(elty,10),C)
+ @test_throws DimensionMismatch LAPACK.ptsv!(rdv,ev,ones(elty,11,11))
  end
 end
 
 @testset "pttrf and pttrs" begin
  @testset for elty in (Float32, Float64, Complex64, Complex128)
- dv = real(ones(elty,10))
+ dv = ones(elty,10)
  ev = zeros(elty,9)
+ rdv = real(dv)
  A = SymTridiagonal(dv,ev)
  if elty <: Complex
  A = Tridiagonal(conj(ev),dv,ev)
  end
- dv,ev = LAPACK.pttrf!(dv,ev)
- @test_throws DimensionMismatch LAPACK.pttrf!(dv,ones(elty,10))
+ rdv,ev = LAPACK.pttrf!(rdv,ev)
+ @test_throws DimensionMismatch LAPACK.pttrf!(rdv,dv)
  B = rand(elty,10,10)
  C = copy(B)
  if elty <: Complex
- @test A\B ≈ LAPACK.pttrs!('U',dv,ev,C)
- @test_throws DimensionMismatch LAPACK.pttrs!('U',dv,ones(elty,10),C)
- @test_throws DimensionMismatch LAPACK.pttrs!('U',dv,ev,rand(elty,11,11))
+ @test A\B ≈ LAPACK.pttrs!('U',rdv,ev,C)
+ @test_throws DimensionMismatch LAPACK.pttrs!('U',rdv,ones(elty,10),C)
+ @test_throws DimensionMismatch LAPACK.pttrs!('U',rdv,ev,rand(elty,11,11))
  else
- @test A\B ≈ LAPACK.pttrs!(dv,ev,C)
- @test_throws DimensionMismatch LAPACK.pttrs!(dv,ones(elty,10),C)
- @test_throws DimensionMismatch LAPACK.pttrs!(dv,ev,rand(elty,11,11))
+ @test A\B ≈ LAPACK.pttrs!(rdv,ev,C)
+ @test_throws DimensionMismatch LAPACK.pttrs!(rdv,ones(elty,10),C)
+ @test_throws DimensionMismatch LAPACK.pttrs!(rdv,ev,rand(elty,11,11))
  end
  end
 end

test/linalg/tridiag.jl

@@ -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
@@ -119,7 +132,8 @@ B = randn(n,2)
  @test_throws DimensionMismatch Tldlt\rand(elty,n+1)
  @test size(Tldlt) == size(Ts)
  if elty <: AbstractFloat
- @test typeof(convert(Base.LinAlg.LDLt{Float32},Tldlt)) == Base.LinAlg.LDLt{Float32,SymTridiagonal{elty}}
+ @test typeof(convert(Base.LinAlg.LDLt{Float32},Tldlt)) ==
+ Base.LinAlg.LDLt{Float32,SymTridiagonal{elty,Vector{elty}}}
  end
  for vv in (vs, view(vs, 1:n))
  invFsv = Fs\vv
@@ -216,7 +230,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)

test/show.jl

@@ -597,7 +597,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"