parametrize SymTridiagonal on the wrapped vector type (#23035)

* parametrize SymTridiagonal on the wrapped vector type

Author: fredrikekre

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"