diff --git a/src/FillArrays.jl b/src/FillArrays.jl index 7fddb28e..b9694d8e 100644 --- a/src/FillArrays.jl +++ b/src/FillArrays.jl @@ -11,7 +11,8 @@ import Base: size, getindex, setindex!, IndexStyle, checkbounds, convert, import LinearAlgebra: rank, svdvals!, tril, triu, tril!, triu!, diag, transpose, adjoint, fill!, dot, norm2, norm1, normInf, normMinusInf, normp, lmul!, rmul!, diagzero, AdjointAbsVec, TransposeAbsVec, - issymmetric, ishermitian, AdjOrTransAbsVec, checksquare, mul!, kron, AbstractTriangular + issymmetric, ishermitian, AdjOrTransAbsVec, checksquare, mul!, kron, AbstractTriangular, + eigvecs, eigvals, eigen import Base.Broadcast: broadcasted, DefaultArrayStyle, broadcast_shape, BroadcastStyle, Broadcasted diff --git a/src/fillalgebra.jl b/src/fillalgebra.jl index f98ae605..838a662c 100644 --- a/src/fillalgebra.jl +++ b/src/fillalgebra.jl @@ -564,3 +564,41 @@ end triu(A::AbstractZerosMatrix, k::Integer=0) = A tril(A::AbstractZerosMatrix, k::Integer=0) = A + +# eigen +_eigenind(λ0, n, sortby) = (isnothing(sortby) || sortby(λ0) <= sortby(zero(λ0))) ? 1 : n + +function eigvals(A::AbstractFillMatrix{<:Union{Real, Complex}}; sortby = nothing) + Base.require_one_based_indexing(A) + n = checksquare(A) + # only one non-trivial eigenvalue for a rank-1 matrix + λ0 = float(getindex_value(A)) * n + ind = _eigenind(λ0, n, sortby) + OneElement(λ0, ind, n) +end + +function eigvecs(A::AbstractFillMatrix{<:Union{Real, Complex}}; sortby = nothing) + Base.require_one_based_indexing(A) + n = checksquare(A) + M = similar(A, real(float(eltype(A)))) + n == 0 && return M + val = getindex_value(A) + ind = _eigenind(val, n, sortby) + # The non-trivial eigenvector is normalize(ones(n)) + M[:, ind] .= inv(sqrt(n)) + # eigenvectors corresponding to zero eigenvalues + for (i, j) in enumerate(axes(M,2)[(ind == 1) .+ (1:end-1)]) + # The eigenvectors are v = normalize([ones(n-1); -(n-1)]), and sum(v) == 0 + # The ordering is arbitrary, + # and we choose to order in terms of the number of non-zero elements + nrm = 1/sqrt(i*(i+1)) + M[1:i, j] .= nrm + M[i+1, j] = -i * nrm + M[i+2:end, j] .= zero(eltype(M)) + end + return M +end + +function eigen(A::AbstractFillMatrix{<:Union{Real, Complex}}; sortby = nothing) + Eigen(eigvals(A; sortby), eigvecs(A; sortby)) +end diff --git a/test/runtests.jl b/test/runtests.jl index b3474f05..6662a9b2 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -2984,6 +2984,32 @@ end @test tril(Z, 2) === Z end +@testset "eigen" begin + sortby = x -> (real(x), imag(x)) + @testset "AbstractFill" begin + sizes = VERSION >= v"1.10" ? (0, 1, 4) : (1, 4) + @testset for val in (2.0, -2, 3+2im, 4 - 5im, 2im), n in sizes + sortby_val = iszero(real(val)) ? imag : sortby + F = Fill(val, n, n) + M = Matrix(F) + @test eigvals(F; sortby = sortby_val) ≈ eigvals(M; sortby = sortby_val) + λ, V = eigen(F; sortby = sortby_val) + @test λ == eigvals(F; sortby = sortby_val) + @test V'V ≈ I + @test F * V ≈ V * Diagonal(λ) + end + @testset for MT in (Ones, Zeros), T in (Float64, Int, ComplexF64), n in sizes + F = MT{T}(n,n) + M = Matrix(F) + @test eigvals(F; sortby) ≈ eigvals(M; sortby) + λ, V = eigen(F; sortby) + @test λ == eigvals(F; sortby) + @test V'V ≈ I + @test F * V ≈ V * Diagonal(λ) + end + end +end + @testset "Diagonal conversion (#389)" begin @test convert(Diagonal{Int, Vector{Int}}, Zeros(5,5)) isa Diagonal{Int,Vector{Int}} @test convert(Diagonal{Int, Vector{Int}}, Zeros(5,5)) == zeros(5,5)