parametrize SymTridiagonal on the wrapped vector type

Author: fredrikekre

NEWS.md

@@ -73,8 +73,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` type definition has changed from `Diagonal{T}` to
-    `Diagonal{T,V<:AbstractVector{T}}` ([#22718]).
+  * 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]).
 
@@ -142,8 +144,9 @@ Library improvements
 
   * `Char`s can now be concatenated with `String`s and/or other `Char`s using `*` ([#22532]).
 
-  * `Diagonal` is now parameterized on the type of the wrapped vector. This allows
-    for `Diagonal` matrices with arbitrary `AbstractVector`s ([#22718]).
+  * `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]).

base/linalg/tridiag.jl

@@ -3,55 +3,80 @@
 #### 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)
 
-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{T}, ev::AbstractVector{S}) where {T,S}
+    R = promote_type(T, S)
+    SymTridiagonal(convert(AbstractVector{R}, dv), convert(AbstractVector{R}, ev))
+end
 
-function SymTridiagonal(dv::AbstractVector{Td}, ev::AbstractVector{Te}) where {Td,Te}
-    T = promote_type(Td,Te)
+function SymTridiagonal(dv::AbstractVector{T}, ev::AbstractVector{T}) where T
     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 `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}:
+ 1  2  ⋅
+ 2  4  5
+ ⋅  5  6
+```
+"""
 function SymTridiagonal(A::AbstractMatrix)
     if diag(A,1) == diag(A,-1)
         SymTridiagonal(diag(A), diag(A,1))
@@ -61,10 +86,10 @@ function SymTridiagonal(A::AbstractMatrix)
 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 +336,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 +366,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 +660,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

base/test.jl

@@ -1315,6 +1315,10 @@ end
 GenericArray{T}(args...) where {T} = GenericArray(Array{T}(args...))
 GenericArray{T,N}(args...) where {T,N} = GenericArray(Array{T,N}(args...))
 
+function Base.convert(::Type{AbstractArray{T,N}}, a::GenericArray) where {T,N}
+    GenericArray{T,N}(convert(Array{T,N}, a.a))
+end
+
 Base.eachindex(a::GenericArray) = eachindex(a.a)
 Base.indices(a::GenericArray) = indices(a.a)
 Base.length(a::GenericArray) = length(a.a)

test/linalg/tridiag.jl

@@ -28,6 +28,15 @@ 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)
+            @test ST::SymTridiagonal{elty, typeof(x)} == Matrix(ST)
+            @test ST.dv === x
+            @test ST.ev === y
+        end
+    end
     ε = eps(abs2(float(one(elty))))
     T = Tridiagonal(dl, d, du)
     Ts = SymTridiagonal(d, dl)
@@ -225,6 +234,26 @@ for elty in (Float32, Float64, Complex64, Complex128, Int)
     end
 end
 
+@testset "arbitrary AbstractVector in SymTridiagonal" begin
+    S, T = Float32, Float64
+    R = promote_type(T, S)
+    dvT = rand(T, n);   dvS = rand(S, n)
+    evT = rand(T, n-1); evS = rand(S, n-1)
+    GdvT = GenericArray(dvT); GdvS = GenericArray(dvS)
+    GevT = GenericArray(evT); GevS = GenericArray(evS)
+    @test isa(SymTridiagonal(dvT, evT), SymTridiagonal{T,Vector{T}})
+    @test isa(SymTridiagonal(dvT, evS), SymTridiagonal{R,Vector{R}})
+    @test isa(SymTridiagonal(dvS, evT), SymTridiagonal{R,Vector{R}})
+
+    @test isa(SymTridiagonal(GdvT, GevT), SymTridiagonal{T,GenericArray{T,1}})
+    @test isa(SymTridiagonal(GdvT, GevS), SymTridiagonal{R,GenericArray{R,1}})
+    @test isa(SymTridiagonal(GdvS, GevT), SymTridiagonal{R,GenericArray{R,1}})
+
+    @test isa(SymTridiagonal(GdvT, evT), SymTridiagonal{T,Vector{T}})
+    @test isa(SymTridiagonal(GdvT, evS), SymTridiagonal{R,Vector{R}})
+    @test isa(SymTridiagonal(GdvS, evT), SymTridiagonal{R,Vector{R}})
+end
+
 #Test equivalence of eigenvectors/singular vectors taking into account possible phase (sign) differences
 function test_approx_eq_vecs(a::StridedVecOrMat{S}, b::StridedVecOrMat{T}, error=nothing) where {S<:Real,T<:Real}
     n = size(a, 1)

test/show.jl

@@ -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}:\n 1  5   ⋅   ⋅\n ⋅  6  10   ⋅\n ⋅  ⋅  11  15\n ⋅  ⋅   ⋅  16"
 @test replstr(Bidiagonal(A,:L)) == "4×4 Bidiagonal{$Int}:\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"