complex skew tridiagonal

[stevengj]

Your current SkewHermTridiagonal only stores the lower diagonal ev, but this is only correct for real skew-tridiagonal. For complex skew-tridiagonal matrices, the diagonal may be nonzero (albeit purely imaginary).

To fix this, you need to add an optional field to store the diagonal imaginary parts if T is complex. Something like:

struct SkewHermTridiagonal{T, V<:AbstractVector{T}, Vim<:Union{AbstractVector{<:Real},Nothing}} <: AbstractMatrix{T}
    ev::V                        # subdiagonal
    dvim::Vim               # diagonal imaginary parts (may be nothing if T is real)
    function SkewHermTridiagonal{T, V, Vim}(ev, dvim) where {T, V<:AbstractVector{T}, Vim}
        LA.require_one_based_indexing(ev)
        if Vim !== Nothing
            LA.require_one_based_indexing(dvim)
            eltype(dvim) === real(T) || throw(ArgumentError("mismatch between $(real(T)) and $(eltype(dvim))"))
        end
        new{T, V, Vim}(ev, dvim)
    end
end

# real skew-symmetric case
SkewHermTridiagonal(ev::AbstractVector{T}) where {T<:Real} = SkewHermTridiagonal{T, typeof(ev), Nothing}(ev, nothing)
SkewHermTridiagonal{T}(ev::AbstractVector) where {T<:Real} =
    SkewHermTridiagonal(convert(AbstractVector{T}, ev)::AbstractVector{T})

# complex skew-hermitian case
SkewHermTridiagonal(ev::AbstractVector{Complex{T}}, dvim::AbstractVector{T}) where T = SkewHermTridiagonal{Complex{T}, typeof(ev), typeof(dvim)}(ev, dvim)
SkewHermTridiagonal{Complex{T}}(ev::AbstractVector, dvim::AbstractVector) where T =
    SkewHermTridiagonal(convert(AbstractVector{Complex{T}}, ev)::AbstractVector{Complex{T}}, convert(AbstractVector{T}, dvim)::AbstractVector{T})

and then you would need to specialize getindex and other functions, e.g.

Base.@propagate_inbounds function Base.getindex(A::SkewHermTridiagonal{T}, i::Integer, j::Integer) where T
    @boundscheck checkbounds(A, i, j)
    if i == j + 1
        return @inbounds A.ev[j]
    elseif i + 1 == j
        return @inbounds -A.ev[i]'
    elseif T <: Complex && i == j
        return complex(zero(real(T)), A.dvim[i])
    else
        return zero(T)
    end
end

(Note the compiler will optimize out the T <: Complex check.)

[smataigne]

I forgot about the imaginary diagonal elements! I will correct this quickly. Thanks a lot!

[smataigne]

This is fixed. I must just correct the copyeigtype function